/[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 1660 by ksteube, Mon Jul 21 03:23: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, 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              int mpi_error;
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    /* optimize node labeling*/      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], mpi_error;
415            temp1[0] = (int) 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 = (ElementTypeId) 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        if (numEle % mpi_info->size != 0) chunkSize++; /* Remainder from numEle / mpi_info->size will be spread out one-per-CPU */
440        if (mpi_info->rank == 0) {  /* Master */
441          for (;;) {            /* Infinite loop */
442            chunkEle = 0;
443            for (i0=0; i0<numEle*(2+numNodes)+1; i0++) tempInts[i0] = -1;
444            for (i0=0; i0<chunkSize; i0++) {
445              if (totalEle >= numEle) break; /* End infinite loop */
446              fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
447              for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
448              fscanf(fileHandle_p, "\n");
449              totalEle++;
450              chunkEle++;
451            }
452            /* Eventually we'll send chunk of nodes to each CPU except 0 itself, here goes one of them */
453            if (nextCPU < mpi_info->size) {
454    #ifdef PASO_MPI
455                  int mpi_error;
456    
457    if (Finley_noError()) {            tempInts[numEle*(2+numNodes)] = chunkEle;
458        if (optimize_labeling) Finley_Mesh_optimizeNodeLabeling(mesh_p);  printf("ksteube CPU=%d/%d send to %d\n", mpi_info->rank, mpi_info->size, nextCPU);
459    }            mpi_error = MPI_Send(tempInts, numEle*(2+numNodes)+1, MPI_INT, nextCPU, 81722, mpi_info->comm);
460              if ( mpi_error != MPI_SUCCESS ) {
461                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Elements failed");
462                    return NULL;
463              }
464    #endif
465              nextCPU++;
466            }
467            if (totalEle >= numEle) break; /* End infinite loop */
468          } /* Infinite loop */
469        }   /* End master */
470        else {  /* Worker */
471    #ifdef PASO_MPI
472          /* Each worker receives two messages */
473          MPI_Status status;
474          int mpi_error;
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      /* Change the typeid from NoType later */
550    
551      mesh_p->FaceElements=Finley_ElementFile_alloc(NoType, mesh_p->order, mesh_p->reduced_order, mpi_info);
552      Finley_ElementFile_allocTable(mesh_p->FaceElements, 0);
553    
554      mesh_p->ContactElements=Finley_ElementFile_alloc(NoType, mesh_p->order, mesh_p->reduced_order, mpi_info);
555      Finley_ElementFile_allocTable(mesh_p->ContactElements, 0);
556    
557      mesh_p->Points=Finley_ElementFile_alloc(NoType, mesh_p->order, mesh_p->reduced_order, mpi_info);
558      Finley_ElementFile_allocTable(mesh_p->Points, 0);
559    
560    #endif
561    
562    #if 0 /* comment out the rest of the un-implemented crap for now */
563            /* get the face elements */
564            if (Finley_noError()) {
565                 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
566                 faceTypeID=Finley_RefElement_getTypeId(element_type);
567                 if (faceTypeID==NoType) {
568                   sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s for face elements",element_type);
569                   Finley_setError(VALUE_ERROR,error_msg);
570                 } else {
571                    mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
572                    if (Finley_noError()) {
573                       Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle);
574                       if (Finley_noError()) {
575                          mesh_p->FaceElements->minColor=0;
576                          mesh_p->FaceElements->maxColor=numEle-1;
577                          for (i0 = 0; i0 < numEle; i0++) {
578                            fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);
579                            mesh_p->FaceElements->Color[i0]=i0;
580                            for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {
581                                 fscanf(fileHandle_p, " %d",
582                                    &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);
583                            }   /* for i1 */
584                            fscanf(fileHandle_p, "\n");
585                          } /* for i0 */
586                       }
587                    }
588                 }
589            }
590            /* get the Contact face element */
591            if (Finley_noError()) {
592                 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
593                 contactTypeID=Finley_RefElement_getTypeId(element_type);
594                 if (contactTypeID==NoType) {
595                   sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type);
596                   Finley_setError(VALUE_ERROR,error_msg);
597                 } else {
598                   mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
599                   if (Finley_noError()) {
600                       Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);
601                       if (Finley_noError()) {
602                          mesh_p->ContactElements->minColor=0;
603                          mesh_p->ContactElements->maxColor=numEle-1;
604                          for (i0 = 0; i0 < numEle; i0++) {
605                            fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]);
606                            mesh_p->ContactElements->Color[i0]=i0;
607                            for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {
608                                fscanf(fileHandle_p, " %d",
609                                   &mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);
610                            }   /* for i1 */
611                            fscanf(fileHandle_p, "\n");
612                          } /* for i0 */
613                      }
614                   }
615                 }
616            }
617            /* get the nodal element */
618            if (Finley_noError()) {
619                 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
620                 pointTypeID=Finley_RefElement_getTypeId(element_type);
621                 if (pointTypeID==NoType) {
622                   sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for points",element_type);
623                   Finley_setError(VALUE_ERROR,error_msg);
624                 }
625                 mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
626                 if (Finley_noError()) {
627                    Finley_ElementFile_allocTable(mesh_p->Points, numEle);
628                    if (Finley_noError()) {
629                       mesh_p->Points->minColor=0;
630                       mesh_p->Points->maxColor=numEle-1;
631                       for (i0 = 0; i0 < numEle; i0++) {
632                         fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);
633                         mesh_p->Points->Color[i0]=i0;
634                         for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {
635                             fscanf(fileHandle_p, " %d",
636                                &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);
637                         }  /* for i1 */
638                         fscanf(fileHandle_p, "\n");
639                       } /* for i0 */
640                    }
641                 }
642            }
643            /* get the name tags */
644            if (Finley_noError()) {
645               if (feof(fileHandle_p) == 0) {
646                  fscanf(fileHandle_p, "%s\n", name);
647                  while (feof(fileHandle_p) == 0) {
648                       fscanf(fileHandle_p, "%s %d\n", name, &tag_key);
649                       Finley_Mesh_addTagMap(mesh_p,name,tag_key);
650                  }
651               }
652            }
653    #endif /* comment out the rest of the un-implemented crap for now */
654         }
655         /* close file */
656         if (mpi_info->rank == 0) fclose(fileHandle_p);
657    
658         /*   resolve id's : */
659         /* rearrange elements: */
660    
661         if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
662         if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
663         return mesh_p; /* ksteube temp return for debugging */
664    
665         /* that's it */
666         #ifdef Finley_TRACE
667         printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);
668         #endif
669    
670    /* that's it */    /* that's it */
671    #ifdef Finley_TRACE    if (! Finley_noError()) {
672    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);  
673    }    }
674      Paso_MPIInfo_free( mpi_info );
675    return mesh_p;    return mesh_p;
676  }  }

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

  ViewVC Help
Powered by ViewVC 1.1.26