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

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

  ViewVC Help
Powered by ViewVC 1.1.26