/[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 616 by elspeth, Wed Mar 22 02:46:56 2006 UTC revision 1628 by phornby, Fri Jul 11 13:12:46 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;
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];
   Finley_NodeFile *nodes_p=NULL;  
38    double time0=Finley_timer();    double time0=Finley_timer();
39      FILE *fileHandle_p = NULL;
40      ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;
41      
42    Finley_resetError();    Finley_resetError();
43    
44    /* get file handle */    if (mpi_info->size > 1) {
45    FILE * 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,"%s: Opening file %s for reading failed.",__FILE__,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);  
   nodes_p=Finley_NodeFile_alloc(numDim);  
   if (! Finley_noError()) return NULL;  
   Finley_NodeFile_allocTable(nodes_p, numNodes);  
   if (! Finley_noError()) return NULL;  
   
   if (1 == numDim) {  
       for (i0 = 0; i0 < numNodes; i0++)  
     fscanf(fileHandle_p, "%d %d %d %le\n", &nodes_p->Id[i0],  
            &nodes_p->degreeOfFreedom[i0], &nodes_p->Tag[i0],  
            &nodes_p->Coordinates[INDEX2(0,i0,numDim)]);  
   } else if (2 == numDim) {  
       for (i0 = 0; i0 < numNodes; i0++)  
     fscanf(fileHandle_p, "%d %d %d %le %le\n", &nodes_p->Id[i0],  
            &nodes_p->degreeOfFreedom[i0], &nodes_p->Tag[i0],  
            &nodes_p->Coordinates[INDEX2(0,i0,numDim)],  
            &nodes_p->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", &nodes_p->Id[i0],  
            &nodes_p->degreeOfFreedom[i0], &nodes_p->Tag[i0],  
            &nodes_p->Coordinates[INDEX2(0,i0,numDim)],  
            &nodes_p->Coordinates[INDEX2(1,i0,numDim)],  
            &nodes_p->Coordinates[INDEX2(2,i0,numDim)]);  
   } /* if else else */  
   
   /* get the element type */  
   
   fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);  
   ElementTypeId typeID=Finley_RefElement_getTypeId(element_type);  
   if (typeID==NoType) {  
     sprintf(error_msg,"%s :Unidentified element type %s",__FILE__,element_type);  
     Finley_setError(VALUE_ERROR,error_msg);  
     return NULL;  
228    }    }
229      Paso_MPIInfo_free( mpi_info );
230      return mesh_p;
231    }
232    
233    /* allocate mesh */  Finley_Mesh* Finley_Mesh_read_MPI(char* fname,index_t order, index_t reduced_order,  bool_t optimize)
   
   Finley_Mesh * mesh_p =Finley_Mesh_alloc(name,numDim,order);  
   if (! Finley_noError()) return NULL;  
   mesh_p->Nodes=nodes_p;  
   
   /* read the elements */  
   mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order);  
   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);  
   ElementTypeId faceTypeID=Finley_RefElement_getTypeId(element_type);  
   if (faceTypeID==NoType) {  
     sprintf(error_msg,"%s :Unidentified element type %s for face elements",__FILE__,element_type);  
     Finley_setError(VALUE_ERROR,error_msg);  
     return NULL;  
   }  
   mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order);  
   Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle);  
   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 */  
   fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);  
   ElementTypeId contactTypeID=Finley_RefElement_getTypeId(element_type);  
   if (contactTypeID==NoType) {  
     sprintf(error_msg,"%s: Unidentified element type %s for contact elements",__FILE__,element_type);  
     Finley_setError(VALUE_ERROR,error_msg);  
     return NULL;  
   }  
   mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order);  
   Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);  
   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 */  
   
   /* get the nodal element */  
   fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);  
   ElementTypeId pointTypeID=Finley_RefElement_getTypeId(element_type);  
   if (pointTypeID==NoType) {  
     sprintf(error_msg,"%s: Unidentified element type %s for points",__FILE__,element_type);  
     Finley_setError(VALUE_ERROR,error_msg);  
     return NULL;  
   }  
   mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order);  
234    
235    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 */  
236    
237      Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
238      dim_t numNodes, numDim, numEle, i0, i1;
239      Finley_Mesh *mesh_p=NULL;
240      char name[LenString_MAX],element_type[LenString_MAX],frm[20];
241      char error_msg[LenErrorMsg_MAX];
242      double time0=Finley_timer();
243      FILE *fileHandle_p = NULL;
244      ElementTypeId typeID;
245    
246    /* close file */    /* these are in unimplemented code below */
247    #if 0
248      index_t tag_key;
249      ElementTypeId faceTypeID, contactTypeID, pointTypeID;
250    #endif
251    
252    fclose(fileHandle_p);    Finley_resetError();
253    
254    /*   resolve id's : */    if (mpi_info->rank == 0) {
255         /* get file handle */
256         fileHandle_p = fopen(fname, "r");
257         if (fileHandle_p==NULL) {
258           sprintf(error_msg,"Finley_Mesh_read: Opening file %s for reading failed.",fname);
259           Finley_setError(IO_ERROR,error_msg);
260           Paso_MPIInfo_free( mpi_info );
261           return NULL;
262         }
263    
264         /* read header */
265         sprintf(frm,"%%%d[^\n]",LenString_MAX-1);
266         fscanf(fileHandle_p, frm, name);
267    
268    Finley_Mesh_resolveNodeIds(mesh_p);       /* get the number of nodes */
269         fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);
270      }
271    
272    /* rearrange elements: */  #ifdef PASO_MPI
273      /* MPI Broadcast numDim, numNodes, name */
274      if (mpi_info->size > 0) {
275        int temp1[3], error_code;
276        temp1[0] = numDim;
277        temp1[1] = numNodes;
278        temp1[2] = strlen(name) + 1;
279        error_code = MPI_Bcast (temp1, 3, MPI_INT,  0, mpi_info->comm);
280        if (error_code != MPI_SUCCESS) {
281          Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of temp1 failed");
282          return NULL;
283        }
284        numDim = temp1[0];
285        numNodes = temp1[1];
286        error_code = MPI_Bcast (name, temp1[2], MPI_CHAR, 0, mpi_info->comm);
287        if (error_code != MPI_SUCCESS) {
288          Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of name failed");
289          return NULL;
290        }
291      }
292    #endif
293    
294    Finley_Mesh_prepare(mesh_p);       /* allocate mesh */
295         mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
296         if (Finley_noError()) {
297        int chunkSize = numNodes / mpi_info->size + 1, totalNodes=0, chunkNodes=0, chunkEle=0, nextCPU=1;
298    #ifdef PASO_MPI
299        int mpi_error;
300    #endif
301        int *tempInts = TMPMEMALLOC(numNodes*3+1, index_t);
302        double *tempCoords = TMPMEMALLOC(numNodes*numDim, double);
303    
304        /*
305          Read a chunk of nodes, send to worker CPU if available, copy chunk into local mesh_p
306          It doesn't matter that a CPU has the wrong nodes for its elements, this is sorted out later
307          First chunk sent to CPU 1, second to CPU 2, ...
308          Last chunk stays on CPU 0 (the master)
309          The three columns of integers (Id, gDOF, Tag) are gathered into a single array tempInts and sent together in a single MPI message
310        */
311    
312        if (mpi_info->rank == 0) {  /* Master */
313          for (;;) {            /* Infinite loop */
314            chunkNodes = 0;
315            for (i0=0; i0<numNodes*3+1; i0++) tempInts[i0] = -1;
316            for (i0=0; i0<numNodes*numDim; i0++) tempCoords[i0] = -1.0;
317            for (i1=0; i1<chunkSize; i1++) {
318              if (totalNodes >= numNodes) break;
319                  if (1 == numDim)
320            fscanf(fileHandle_p, "%d %d %d %le\n",
321              &tempInts[0+i1], &tempInts[numNodes+i1], &tempInts[numNodes*2+i1],
322              &tempCoords[i1*numDim+0]);
323                  if (2 == numDim)
324            fscanf(fileHandle_p, "%d %d %d %le %le\n",
325              &tempInts[0+i1], &tempInts[numNodes+i1], &tempInts[numNodes*2+i1],
326              &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1]);
327                  if (3 == numDim)
328            fscanf(fileHandle_p, "%d %d %d %le %le %le\n",
329              &tempInts[0+i1], &tempInts[numNodes+i1], &tempInts[numNodes*2+i1],
330              &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1], &tempCoords[i1*numDim+2]);
331              totalNodes++;
332              chunkNodes++;
333            }
334            /* Eventually we'll send chunk of nodes to each CPU numbered 1 ... mpi_info->size-1, here goes one of them */
335            if (nextCPU < mpi_info->size) {
336    #ifdef PASO_MPI
337              tempInts[numNodes*3] = chunkNodes;
338              /* ksteube The size of this message can and should be brought down to chunkNodes*3+1, must re-org tempInts */
339              mpi_error = MPI_Send(tempInts, numNodes*3+1, MPI_INT, nextCPU, 81720, mpi_info->comm);
340              if ( mpi_error != MPI_SUCCESS ) {
341                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts failed");
342                    return NULL;
343              }
344              mpi_error = MPI_Send(tempCoords, numNodes*numDim, MPI_DOUBLE, nextCPU, 81721, mpi_info->comm);
345              if ( mpi_error != MPI_SUCCESS ) {
346                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempCoords failed");
347                    return NULL;
348              }
349    #endif
350              nextCPU++;
351            }
352            if (totalNodes >= numNodes) break;
353          } /* Infinite loop */
354        }   /* End master */
355        else {  /* Worker */
356    #ifdef PASO_MPI
357          /* Each worker receives two messages */
358          MPI_Status status;
359          mpi_error = MPI_Recv(tempInts, numNodes*3+1, MPI_INT, 0, 81720, mpi_info->comm, &status);
360          if ( mpi_error != MPI_SUCCESS ) {
361                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts failed");
362                return NULL;
363          }
364          mpi_error = MPI_Recv(tempCoords, numNodes*numDim, MPI_DOUBLE, 0, 81721, mpi_info->comm, &status);
365          if ( mpi_error != MPI_SUCCESS ) {
366                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempCoords failed");
367                return NULL;
368          }
369          chunkNodes = tempInts[numNodes*3];
370    #endif
371        }   /* Worker */
372    
373    #if 0
374            /* Display the temp mem for debugging */
375            printf("ksteube tempInts totalNodes=%d:\n", totalNodes);
376            for (i0=0; i0<numNodes*3; i0++) {
377              printf(" %2d", tempInts[i0]);
378              if (i0%numNodes==numNodes-1) printf("\n");
379            }
380            printf("ksteube tempCoords:\n");
381            for (i0=0; i0<chunkNodes*numDim; i0++) {
382              printf(" %20.15e", tempCoords[i0]);
383              if (i0%numDim==numDim-1) printf("\n");
384            }
385    #endif
386    
387        printf("ksteube chunkNodes=%d numNodes=%d\n", chunkNodes, numNodes);
388        /* Copy node data from tempMem to mesh_p */
389            Finley_NodeFile_allocTable(mesh_p->Nodes, chunkNodes);
390            if (Finley_noError()) {
391          for (i0=0; i0<chunkNodes; i0++) {
392            mesh_p->Nodes->Id[i0]               = tempInts[0+i0];
393            mesh_p->Nodes->globalDegreesOfFreedom[i0]       = tempInts[numNodes+i0];
394            mesh_p->Nodes->Tag[i0]              = tempInts[numNodes*2+i0];
395            for (i1=0; i1<numDim; i1++) {
396              mesh_p->Nodes->Coordinates[INDEX2(i1,i0,numDim)]  = tempCoords[i0*numDim+i1];
397            }
398              }
399            }
400    
401        TMPMEMFREE(tempInts);
402        TMPMEMFREE(tempCoords);
403    
404            /* read elements */
405    
406        /* Read the element typeID */
407            if (Finley_noError()) {
408          if (mpi_info->rank == 0) {
409                fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
410                typeID=Finley_RefElement_getTypeId(element_type);
411          }
412    #ifdef PASO_MPI
413          if (mpi_info->size > 0) {
414            int temp1[3];
415            temp1[0] = typeID;
416            temp1[1] = numEle;
417            mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
418            if (mpi_error != MPI_SUCCESS) {
419              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
420              return NULL;
421            }
422            typeID = temp1[0];
423            numEle = temp1[1];
424          }
425    #endif
426              if (typeID==NoType) {
427                sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
428                Finley_setError(VALUE_ERROR, error_msg);
429              }
430        }
431    
432          /* Read the element data */
433          mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
434          numNodes = mesh_p->Elements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
435    
436          if (Finley_noError()) {
437        int *tempInts = TMPMEMALLOC(numEle*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
438        int chunkSize = numEle / mpi_info->size, totalEle=0, nextCPU=1;
439    #ifdef PASO_MPI
440        int mpi_error;
441    #endif
442        if (numEle % mpi_info->size != 0) chunkSize++; /* Remainder from numEle / mpi_info->size will be spread out one-per-CPU */
443        if (mpi_info->rank == 0) {  /* Master */
444          for (;;) {            /* Infinite loop */
445            chunkEle = 0;
446            for (i0=0; i0<numEle*(2+numNodes)+1; i0++) tempInts[i0] = -1;
447            for (i0=0; i0<chunkSize; i0++) {
448              if (totalEle >= numEle) break; /* End infinite loop */
449              fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
450              for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
451              fscanf(fileHandle_p, "\n");
452              totalEle++;
453              chunkEle++;
454            }
455            /* Eventually we'll send chunk of nodes to each CPU except 0 itself, here goes one of them */
456            if (nextCPU < mpi_info->size) {
457    #ifdef PASO_MPI
458              tempInts[numEle*(2+numNodes)] = chunkEle;
459    printf("ksteube CPU=%d/%d send to %d\n", mpi_info->rank, mpi_info->size, nextCPU);
460              mpi_error = MPI_Send(tempInts, numEle*(2+numNodes)+1, MPI_INT, nextCPU, 81722, mpi_info->comm);
461              if ( mpi_error != MPI_SUCCESS ) {
462                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Elements failed");
463                    return NULL;
464              }
465    #endif
466              nextCPU++;
467            }
468            if (totalEle >= numEle) break; /* End infinite loop */
469          } /* Infinite loop */
470        }   /* End master */
471        else {  /* Worker */
472    #ifdef PASO_MPI
473          /* Each worker receives two messages */
474          MPI_Status status;
475    printf("ksteube CPU=%d/%d recv on %d\n", mpi_info->rank, mpi_info->size, mpi_info->rank);
476          mpi_error = MPI_Recv(tempInts, numEle*(2+numNodes)+1, MPI_INT, 0, 81722, mpi_info->comm, &status);
477          if ( mpi_error != MPI_SUCCESS ) {
478                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for Elements failed");
479                return NULL;
480          }
481          chunkEle = tempInts[numEle*(2+numNodes)];
482    #endif
483        }   /* Worker */
484    #if 1
485        /* Display the temp mem for debugging */
486        printf("ksteube tempInts numEle=%d chunkEle=%d AAA:\n", numEle, chunkEle);
487        for (i0=0; i0<numEle*(numNodes+2); i0++) {
488          printf(" %2d", tempInts[i0]);
489          if (i0%(numNodes+2)==numNodes+2-1) printf("\n");
490        }
491    #endif
492    
493        /* Copy Element data from tempInts to mesh_p */
494        Finley_ElementFile_allocTable(mesh_p->Elements, chunkEle);
495            mesh_p->Elements->minColor=0;
496            mesh_p->Elements->maxColor=chunkEle-1;
497            if (Finley_noError()) {
498              #pragma omp parallel for private (i0, i1)
499          for (i0=0; i0<chunkEle; i0++) {
500            mesh_p->Elements->Id[i0]    = tempInts[i0*(2+numNodes)+0];
501            mesh_p->Elements->Tag[i0]   = tempInts[i0*(2+numNodes)+1];
502            for (i1 = 0; i1 < numNodes; i1++) {
503              mesh_p->Elements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
504            }
505              }
506            }
507    
508        TMPMEMFREE(tempInts);
509          }
510    
511    
512    #if 0 /* this is the original code for reading elements */
513            /* read elements */
514            if (Finley_noError()) {
515    
516               fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
517               typeID=Finley_RefElement_getTypeId(element_type);
518               if (typeID==NoType) {
519                 sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s",element_type);
520                 Finley_setError(VALUE_ERROR,error_msg);
521               } else {
522                 /* read the elements */
523                 mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
524                 if (Finley_noError()) {
525                     Finley_ElementFile_allocTable(mesh_p->Elements, numEle);
526                     mesh_p->Elements->minColor=0;
527                     mesh_p->Elements->maxColor=numEle-1;
528                     if (Finley_noError()) {
529                        for (i0 = 0; i0 < numEle; i0++) {
530                          fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]);
531                          mesh_p->Elements->Color[i0]=i0;
532                          for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {
533                               fscanf(fileHandle_p, " %d",
534                                  &mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);
535                          } /* for i1 */
536                          fscanf(fileHandle_p, "\n");
537                        } /* for i0 */
538                     }
539                 }
540              }
541            }
542    #endif
543    
544    printf("ksteube CPU=%d/%d Element typeID=%d\n", mpi_info->rank, mpi_info->size, typeID);
545    
546    #if 1
547    
548      /* Define other structures to keep mesh_write from crashing */
549    
550      mesh_p->FaceElements=Finley_ElementFile_alloc(0, mesh_p->order, mesh_p->reduced_order, mpi_info);
551      Finley_ElementFile_allocTable(mesh_p->FaceElements, 0);
552    
553      mesh_p->ContactElements=Finley_ElementFile_alloc(0, mesh_p->order, mesh_p->reduced_order, mpi_info);
554      Finley_ElementFile_allocTable(mesh_p->ContactElements, 0);
555    
556      mesh_p->Points=Finley_ElementFile_alloc(0, mesh_p->order, mesh_p->reduced_order, mpi_info);
557      Finley_ElementFile_allocTable(mesh_p->Points, 0);
558    
559    #endif
560    
561    #if 0 /* comment out the rest of the un-implemented crap for now */
562            /* get the face elements */
563            if (Finley_noError()) {
564                 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
565                 faceTypeID=Finley_RefElement_getTypeId(element_type);
566                 if (faceTypeID==NoType) {
567                   sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s for face elements",element_type);
568                   Finley_setError(VALUE_ERROR,error_msg);
569                 } else {
570                    mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
571                    if (Finley_noError()) {
572                       Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle);
573                       if (Finley_noError()) {
574                          mesh_p->FaceElements->minColor=0;
575                          mesh_p->FaceElements->maxColor=numEle-1;
576                          for (i0 = 0; i0 < numEle; i0++) {
577                            fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);
578                            mesh_p->FaceElements->Color[i0]=i0;
579                            for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {
580                                 fscanf(fileHandle_p, " %d",
581                                    &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);
582                            }   /* for i1 */
583                            fscanf(fileHandle_p, "\n");
584                          } /* for i0 */
585                       }
586                    }
587                 }
588            }
589            /* get the Contact face element */
590            if (Finley_noError()) {
591                 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
592                 contactTypeID=Finley_RefElement_getTypeId(element_type);
593                 if (contactTypeID==NoType) {
594                   sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type);
595                   Finley_setError(VALUE_ERROR,error_msg);
596                 } else {
597                   mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
598                   if (Finley_noError()) {
599                       Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);
600                       if (Finley_noError()) {
601                          mesh_p->ContactElements->minColor=0;
602                          mesh_p->ContactElements->maxColor=numEle-1;
603                          for (i0 = 0; i0 < numEle; i0++) {
604                            fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]);
605                            mesh_p->ContactElements->Color[i0]=i0;
606                            for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {
607                                fscanf(fileHandle_p, " %d",
608                                   &mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);
609                            }   /* for i1 */
610                            fscanf(fileHandle_p, "\n");
611                          } /* for i0 */
612                      }
613                   }
614                 }
615            }
616            /* get the nodal element */
617            if (Finley_noError()) {
618                 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
619                 pointTypeID=Finley_RefElement_getTypeId(element_type);
620                 if (pointTypeID==NoType) {
621                   sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for points",element_type);
622                   Finley_setError(VALUE_ERROR,error_msg);
623                 }
624                 mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
625                 if (Finley_noError()) {
626                    Finley_ElementFile_allocTable(mesh_p->Points, numEle);
627                    if (Finley_noError()) {
628                       mesh_p->Points->minColor=0;
629                       mesh_p->Points->maxColor=numEle-1;
630                       for (i0 = 0; i0 < numEle; i0++) {
631                         fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);
632                         mesh_p->Points->Color[i0]=i0;
633                         for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {
634                             fscanf(fileHandle_p, " %d",
635                                &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);
636                         }  /* for i1 */
637                         fscanf(fileHandle_p, "\n");
638                       } /* for i0 */
639                    }
640                 }
641            }
642            /* get the name tags */
643            if (Finley_noError()) {
644               if (feof(fileHandle_p) == 0) {
645                  fscanf(fileHandle_p, "%s\n", name);
646                  while (feof(fileHandle_p) == 0) {
647                       fscanf(fileHandle_p, "%s %d\n", name, &tag_key);
648                       Finley_Mesh_addTagMap(mesh_p,name,tag_key);
649                  }
650               }
651            }
652    #endif /* comment out the rest of the un-implemented crap for now */
653         }
654         /* close file */
655         if (mpi_info->rank == 0) fclose(fileHandle_p);
656    
657         /*   resolve id's : */
658         /* rearrange elements: */
659    
660         if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
661         if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
662         return mesh_p; /* ksteube temp return for debugging */
663    
664         /* that's it */
665         #ifdef Finley_TRACE
666         printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);
667         #endif
668    
669    /* that's it */    /* that's it */
670    #ifdef Finley_TRACE    if (! Finley_noError()) {
671    printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);         Finley_Mesh_free(mesh_p);
672    #endif    }
673    if (! Finley_noError()) Finley_Mesh_dealloc(mesh_p);    Paso_MPIInfo_free( mpi_info );
674    return mesh_p;    return mesh_p;
675  }  }
 /*  
 * $Log$  
 * Revision 1.5  2005/09/15 03:44:22  jgs  
 * Merge of development branch dev-02 back to main trunk on 2005-09-15  
 *  
 * Revision 1.4  2005/09/01 03:31:36  jgs  
 * Merge of development branch dev-02 back to main trunk on 2005-09-01  
 *  
 * Revision 1.3.2.3  2005/09/08 08:28:39  gross  
 * some cleanup in savevtk  
 *  
 * Revision 1.3.2.2  2005/09/07 06:26:19  gross  
 * the solver from finley are put into the standalone package paso now  
 *  
 * Revision 1.3.2.1  2005/08/24 02:02:18  gross  
 * timing output switched off. solver output can be swiched through getSolution(verbose=True) now.  
 *  
 * Revision 1.3  2005/07/22 03:53:08  jgs  
 * Merge of development branch back to main trunk on 2005-07-22  
 *  
 * Revision 1.2  2005/07/08 04:07:54  jgs  
 * Merge of development branch back to main trunk on 2005-07-08  
 *  
 * Revision 1.1.1.1.2.2  2005/07/18 10:34:54  gross  
 * some informance improvements when reading meshes  
 *  
 * Revision 1.1.1.1.2.1  2005/06/29 02:34:53  gross  
 * some changes towards 64 integers in finley  
 *  
 * Revision 1.1.1.1  2004/10/26 06:53:57  jgs  
 * initial import of project esys2  
 *  
 * Revision 1.1.1.1  2004/06/24 04:00:40  johng  
 * Initial version of eys using boost-python.  
 *  
 *  
 */  
   

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

  ViewVC Help
Powered by ViewVC 1.1.26