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

Legend:
Removed from v.1156  
changed lines
  Added in v.1646

  ViewVC Help
Powered by ViewVC 1.1.26