/[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

trunk/finley/src/Mesh_read.c revision 1028 by gross, Wed Mar 14 00:15:24 2007 UTC temp/finley/src/Mesh_read.c revision 1387 by trankine, Fri Jan 11 07:45:26 2008 UTC
# Line 1  Line 1 
 /*  
  ************************************************************  
  *          Copyright 2006 by ACcESS MNRF                   *  
  *                                                          *  
  *              http://www.access.edu.au                    *  
  *       Primary Business: Queensland, Australia            *  
  *  Licensed under the Open Software License version 3.0    *  
  *     http://www.opensource.org/licenses/osl-3.0.php       *  
  *                                                          *  
  ************************************************************  
 */  
1    
2  /**************************************************************/  /* $Id$ */
3    
4  /*   Finley: read mesh */  /*******************************************************
5     *
6     *           Copyright 2003-2007 by ACceSS MNRF
7     *       Copyright 2007 by University of Queensland
8     *
9     *                http://esscc.uq.edu.au
10     *        Primary Business: Queensland, Australia
11     *  Licensed under the Open Software License version 3.0
12     *     http://www.opensource.org/licenses/osl-3.0.php
13     *
14     *******************************************************/
15    
16  /**************************************************************/  /**************************************************************/
17    
18  /*   Author: gross@access.edu.au */  /*   Finley: read mesh */
 /*   Version: $Id$ */  
19    
20  /**************************************************************/  /**************************************************************/
21    
# Line 27  Line 25 
25    
26  /*  reads a mesh from a Finley file of name fname */  /*  reads a mesh from a Finley file of name fname */
27    
28  Finley_Mesh* Finley_Mesh_read(char* fname,index_t order) {  Finley_Mesh* Finley_Mesh_read(char* fname,index_t order, index_t reduced_order,  bool_t optimize)
29    
30    {
31    
32      Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
33    dim_t numNodes, numDim, numEle, i0, i1;    dim_t numNodes, numDim, numEle, i0, i1;
34      index_t tag_key;
35    Finley_Mesh *mesh_p=NULL;    Finley_Mesh *mesh_p=NULL;
36    char name[LenString_MAX],element_type[LenString_MAX],frm[20];    char name[LenString_MAX],element_type[LenString_MAX],frm[20];
37    char error_msg[LenErrorMsg_MAX];    char error_msg[LenErrorMsg_MAX];
38    double time0=Finley_timer();    double time0=Finley_timer();
39    FILE *fileHandle_p = NULL;    FILE *fileHandle_p = NULL;
40    ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;    ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;
41      
42    Finley_resetError();    Finley_resetError();
43    
44    /* get file handle */    if (mpi_info->size > 1) {
45    fileHandle_p = fopen(fname, "r");       Finley_setError(SYSTEM_ERROR,"Finley_Mesh_read: MPI is not suporrted yet.");
46    if (fileHandle_p==NULL) {    } else {
47      sprintf(error_msg,"Finley_Mesh_read: Opening file %s for reading failed.",fname);       /* get file handle */
48      Finley_setError(IO_ERROR,error_msg);       fileHandle_p = fopen(fname, "r");
49      return NULL;       if (fileHandle_p==NULL) {
50           sprintf(error_msg,"Finley_Mesh_read: Opening file %s for reading failed.",fname);
51           Finley_setError(IO_ERROR,error_msg);
52           Paso_MPIInfo_free( mpi_info );
53           return NULL;
54         }
55      
56         /* read header */
57         sprintf(frm,"%%%d[^\n]",LenString_MAX-1);
58         fscanf(fileHandle_p, frm, name);
59      
60         /* get the nodes */
61      
62         fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);
63         /* allocate mesh */
64         mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
65         if (Finley_noError()) {
66      
67            /* read nodes */
68            Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
69            if (Finley_noError()) {
70               if (1 == numDim) {
71                   for (i0 = 0; i0 < numNodes; i0++)
72                    fscanf(fileHandle_p, "%d %d %d %le\n", &mesh_p->Nodes->Id[i0],
73                           &mesh_p->Nodes->globalDegreesOfFreedom[i0], &mesh_p->Nodes->Tag[i0],
74                           &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)]);
75               } else if (2 == numDim) {
76                        for (i0 = 0; i0 < numNodes; i0++)
77                              fscanf(fileHandle_p, "%d %d %d %le %le\n", &mesh_p->Nodes->Id[i0],
78                                     &mesh_p->Nodes->globalDegreesOfFreedom[i0], &mesh_p->Nodes->Tag[i0],
79                                     &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],
80                                     &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)]);
81               } else if (3 == numDim) {
82                        for (i0 = 0; i0 < numNodes; i0++)
83                              fscanf(fileHandle_p, "%d %d %d %le %le %le\n", &mesh_p->Nodes->Id[i0],
84                                     &mesh_p->Nodes->globalDegreesOfFreedom[i0], &mesh_p->Nodes->Tag[i0],
85                                     &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],
86                                     &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)],
87                                     &mesh_p->Nodes->Coordinates[INDEX2(2,i0,numDim)]);
88               } /* if else else */
89            }
90            /* read elements */
91            if (Finley_noError()) {
92      
93               fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
94               typeID=Finley_RefElement_getTypeId(element_type);
95               if (typeID==NoType) {
96                 sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s",element_type);
97                 Finley_setError(VALUE_ERROR,error_msg);
98               } else {
99                 /* read the elements */
100                 mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
101                 if (Finley_noError()) {
102                     Finley_ElementFile_allocTable(mesh_p->Elements, numEle);
103                     mesh_p->Elements->minColor=0;
104                     mesh_p->Elements->maxColor=numEle-1;
105                     if (Finley_noError()) {
106                        for (i0 = 0; i0 < numEle; i0++) {
107                          fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]);
108                          mesh_p->Elements->Color[i0]=i0;
109                          for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {
110                               fscanf(fileHandle_p, " %d",
111                                  &mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);
112                          } /* for i1 */
113                          fscanf(fileHandle_p, "\n");
114                        } /* for i0 */
115                     }
116                 }
117              }
118            }
119            /* get the face elements */
120            if (Finley_noError()) {
121                 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
122                 faceTypeID=Finley_RefElement_getTypeId(element_type);
123                 if (faceTypeID==NoType) {
124                   sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s for face elements",element_type);
125                   Finley_setError(VALUE_ERROR,error_msg);
126                 } else {
127                    mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
128                    if (Finley_noError()) {
129                       Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle);
130                       if (Finley_noError()) {
131                          mesh_p->FaceElements->minColor=0;
132                          mesh_p->FaceElements->maxColor=numEle-1;
133                          for (i0 = 0; i0 < numEle; i0++) {
134                            fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);
135                            mesh_p->FaceElements->Color[i0]=i0;
136                            for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {
137                                 fscanf(fileHandle_p, " %d",
138                                    &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);
139                            }   /* for i1 */
140                            fscanf(fileHandle_p, "\n");
141                          } /* for i0 */
142                       }
143                    }
144                 }
145            }
146            /* get the Contact face element */
147            if (Finley_noError()) {
148                 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
149                 contactTypeID=Finley_RefElement_getTypeId(element_type);
150                 if (contactTypeID==NoType) {
151                   sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type);
152                   Finley_setError(VALUE_ERROR,error_msg);
153                 } else {
154                   mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
155                   if (Finley_noError()) {
156                       Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);
157                       if (Finley_noError()) {
158                          mesh_p->ContactElements->minColor=0;
159                          mesh_p->ContactElements->maxColor=numEle-1;
160                          for (i0 = 0; i0 < numEle; i0++) {
161                            fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]);
162                            mesh_p->ContactElements->Color[i0]=i0;
163                            for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {
164                                fscanf(fileHandle_p, " %d",
165                                   &mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);
166                            }   /* for i1 */
167                            fscanf(fileHandle_p, "\n");
168                          } /* for i0 */
169                      }
170                   }
171                 }
172            }  
173            /* get the nodal element */
174            if (Finley_noError()) {
175                 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
176                 pointTypeID=Finley_RefElement_getTypeId(element_type);
177                 if (pointTypeID==NoType) {
178                   sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for points",element_type);
179                   Finley_setError(VALUE_ERROR,error_msg);
180                 }
181                 mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
182                 if (Finley_noError()) {
183                    Finley_ElementFile_allocTable(mesh_p->Points, numEle);
184                    if (Finley_noError()) {
185                       mesh_p->Points->minColor=0;
186                       mesh_p->Points->maxColor=numEle-1;
187                       for (i0 = 0; i0 < numEle; i0++) {
188                         fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);
189                         mesh_p->Points->Color[i0]=i0;
190                         for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {
191                             fscanf(fileHandle_p, " %d",
192                                &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);
193                         }  /* for i1 */
194                         fscanf(fileHandle_p, "\n");
195                       } /* for i0 */
196                    }
197                 }
198            }
199            /* get the name tags */
200            if (Finley_noError()) {
201               if (feof(fileHandle_p) == 0) {
202                  fscanf(fileHandle_p, "%s\n", name);
203                  while (feof(fileHandle_p) == 0) {
204                       fscanf(fileHandle_p, "%s %d\n", name, &tag_key);
205                       Finley_Mesh_addTagMap(mesh_p,name,tag_key);
206                  }
207               }
208            }
209         }
210         /* close file */
211         fclose(fileHandle_p);
212      
213         /*   resolve id's : */
214         /* rearrange elements: */
215      
216         if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
217         if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
218      
219         /* that's it */
220         #ifdef Finley_TRACE
221         printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);
222         #endif
223    }    }
224    
225    /* read header */    /* that's it */
226    sprintf(frm,"%%%d[^\n]",LenString_MAX-1);    if (! Finley_noError()) {
227    fscanf(fileHandle_p, frm, name);         Finley_Mesh_free(mesh_p);
   
   /* get the nodes */  
   
   fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);  
   /* allocate mesh */  
 #ifndef PASO_MPI  
   mesh_p = Finley_Mesh_alloc(name,numDim,order);  
   if (! Finley_noError()) return NULL;  
 #else  
   /* TODO */  
 #endif  
   
 #ifndef PASO_MPI  
   Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);  
   if (! Finley_noError()) return NULL;  
 #else  
   /* TODO */  
 #endif  
   
   if (1 == numDim) {  
       for (i0 = 0; i0 < numNodes; i0++)  
     fscanf(fileHandle_p, "%d %d %d %le\n", &mesh_p->Nodes->Id[i0],  
            &mesh_p->Nodes->degreeOfFreedom[i0], &mesh_p->Nodes->Tag[i0],  
            &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)]);  
   } else if (2 == numDim) {  
       for (i0 = 0; i0 < numNodes; i0++)  
     fscanf(fileHandle_p, "%d %d %d %le %le\n", &mesh_p->Nodes->Id[i0],  
            &mesh_p->Nodes->degreeOfFreedom[i0], &mesh_p->Nodes->Tag[i0],  
            &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],  
            &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)]);  
   } else if (3 == numDim) {  
       for (i0 = 0; i0 < numNodes; i0++)  
     fscanf(fileHandle_p, "%d %d %d %le %le %le\n", &mesh_p->Nodes->Id[i0],  
            &mesh_p->Nodes->degreeOfFreedom[i0], &mesh_p->Nodes->Tag[i0],  
            &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],  
            &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)],  
            &mesh_p->Nodes->Coordinates[INDEX2(2,i0,numDim)]);  
   } /* if else else */  
   
   /* get the element type */  
   
   fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);  
   typeID=Finley_RefElement_getTypeId(element_type);  
   if (typeID==NoType) {  
     sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s",element_type);  
     Finley_setError(VALUE_ERROR,error_msg);  
     return NULL;  
   }  
   /* read the elements */  
 #ifndef PASO_MPI  
   mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order);  
 #else  
   /* TODO */  
 #endif  
   Finley_ElementFile_allocTable(mesh_p->Elements, numEle);  
   mesh_p->Elements->minColor=0;  
   mesh_p->Elements->maxColor=numEle-1;  
   for (i0 = 0; i0 < numEle; i0++) {  
     fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]);  
     mesh_p->Elements->Color[i0]=i0;  
     for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {  
          fscanf(fileHandle_p, " %d",  
             &mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);  
     }   /* for i1 */  
     fscanf(fileHandle_p, "\n");  
   } /* for i0 */  
   
   /* get the face elements */  
   fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);  
   faceTypeID=Finley_RefElement_getTypeId(element_type);  
   faceTypeID=Finley_RefElement_getTypeId(element_type);  
   if (faceTypeID==NoType) {  
     sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s for face elements",element_type);  
     Finley_setError(VALUE_ERROR,error_msg);  
     return NULL;  
   }  
 #ifndef PASO_MPI  
   mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order);  
 #else  
   /* TODO */  
 #endif  
   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;  
   }  
 #ifndef PASO_MPI  
   mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order);  
 #else  
   /* TODO */  
 #endif  
   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  #ifndef PASO_MPI    Paso_MPIInfo_free( mpi_info );
230    mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order);    return mesh_p;
231  #else  }
   /* TODO */  
 #endif  
   Finley_ElementFile_allocTable(mesh_p->Points, numEle);  
   mesh_p->Points->minColor=0;  
   mesh_p->Points->maxColor=numEle-1;  
   for (i0 = 0; i0 < numEle; i0++) {  
     fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);  
     mesh_p->Points->Color[i0]=i0;  
     for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {  
         fscanf(fileHandle_p, " %d",  
            &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);  
     }   /* for i1 */  
     fscanf(fileHandle_p, "\n");  
   } /* for i0 */  
232    
233    Finley_Mesh* Finley_Mesh_read_MPI(char* fname,index_t order, index_t reduced_order,  bool_t optimize)
234    
235    /* close file */  {
236    
237    fclose(fileHandle_p);    Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
238      dim_t numNodes, numDim, numEle, i0, i1;
239      index_t tag_key;
240      Finley_Mesh *mesh_p=NULL;
241      char name[LenString_MAX],element_type[LenString_MAX],frm[20];
242      char error_msg[LenErrorMsg_MAX];
243      double time0=Finley_timer();
244      FILE *fileHandle_p = NULL;
245      ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;
246      
247      Finley_resetError();
248    
249    /*   resolve id's : */    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         }
258      
259         /* read header */
260         sprintf(frm,"%%%d[^\n]",LenString_MAX-1);
261         fscanf(fileHandle_p, frm, name);
262      
263         /* get the nodes */
264      
265         fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);
266      }
267    
268    Finley_Mesh_resolveNodeIds(mesh_p);  #ifdef PASO_MPI
269      /* MPI Broadcast numDim, numNodes, name */
270      if (mpi_info->size > 0) {
271        int temp1[3], error_code;
272        temp1[0] = numDim;
273        temp1[1] = numNodes;
274        temp1[2] = strlen(name) + 1;
275        error_code = MPI_Bcast (temp1, 3, MPI_INT,  0, mpi_info->comm);
276        if (error_code != MPI_SUCCESS) {
277          Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of temp1 failed");
278          return NULL;
279        }
280        numDim = temp1[0];
281        numNodes = temp1[1];
282        error_code = MPI_Bcast (name, temp1[2], MPI_CHAR, 0, mpi_info->comm);
283        if (error_code != MPI_SUCCESS) {
284          Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of name failed");
285          return NULL;
286        }
287      }
288    #endif
289      printf("ksteube CPU=%d/%d name='%s' numDim=%d numNodes=%d\n", mpi_info->rank, mpi_info->size, name, numDim, numNodes);
290    
291    /* rearrange elements: */       /* allocate mesh */
292         mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
293         if (Finley_noError()) {
294        int chunkSize = numNodes / mpi_info->size + 1, totalNodes=0, chunkNodes=0, nextCPU=1, mpi_error;
295        int *tempInts = TMPMEMALLOC(numNodes*3+1, index_t);
296        double *tempCoords = TMPMEMALLOC(numNodes*numDim, double);
297    
298        /*
299          Read a chunk of nodes, send to worker CPU if necessary, copy chunk into local mesh_p
300          It doesn't matter that a CPU has the wrong nodes, this is sorted out later
301          First chunk sent to CPU 1, second to CPU 2, ...
302          Last chunk stays on CPU 0 (the master)
303          The three columns of integers (Id, gDOF, Tag) are gathered into a single array tempInts and sent together in a single MPI message
304        */
305    
306        if (mpi_info->rank == 0) {  /* Master */
307          for (;;) {            /* Infinite loop */
308            chunkNodes = 0;
309            for (i0=0; i0<numNodes*3; i0++) tempInts[i0] = -1;
310            for (i0=0; i0<numNodes*numDim; i0++) tempCoords[i0] = -1.0;
311            for (i1=0; i1<chunkSize; i1++) {
312              if (totalNodes >= numNodes) break;
313                  if (1 == numDim)
314            fscanf(fileHandle_p, "%d %d %d %le\n",
315              &tempInts[0+i1], &tempInts[numNodes+i1], &tempInts[numNodes*2+i1],
316              &tempCoords[i1*numDim+0]);
317                  if (2 == numDim)
318            fscanf(fileHandle_p, "%d %d %d %le %le\n",
319              &tempInts[0+i1], &tempInts[numNodes+i1], &tempInts[numNodes*2+i1],
320              &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1]);
321                  if (3 == numDim)
322            fscanf(fileHandle_p, "%d %d %d %le %le %le\n",
323              &tempInts[0+i1], &tempInts[numNodes+i1], &tempInts[numNodes*2+i1],
324              &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1], &tempCoords[i1*numDim+2]);
325              totalNodes++;
326              chunkNodes++;
327            }
328            /* Eventually we'll send chunk of nodes to each CPU numbered 1 ... mpi_info->size-1, here goes one of them */
329            if (nextCPU < mpi_info->size) {
330    #ifdef PASO_MPI
331              tempInts[numNodes*3] = chunkNodes;
332              mpi_error = MPI_Send(tempInts, numNodes*3+1, MPI_INT, nextCPU, 8727, mpi_info->comm);
333              if ( mpi_error != MPI_SUCCESS ) {
334                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts failed");
335                    return NULL;
336              }
337              mpi_error = MPI_Send(tempCoords, numNodes*numDim, MPI_DOUBLE, nextCPU, 8728, mpi_info->comm);
338              if ( mpi_error != MPI_SUCCESS ) {
339                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempCoords failed");
340                    return NULL;
341              }
342    #endif
343              nextCPU++;
344            }
345            if (totalNodes >= numNodes) break;
346          } /* Infinite loop */
347        }   /* End master */
348        else {  /* Worker */
349    #ifdef PASO_MPI
350          /* Each worker receives two messages */
351          MPI_Status status;
352          mpi_error = MPI_Recv(tempInts, numNodes*3+1, MPI_INT, 0, 8727, mpi_info->comm, &status);
353          if ( mpi_error != MPI_SUCCESS ) {
354                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts failed");
355                return NULL;
356          }
357          mpi_error = MPI_Recv(tempCoords, numNodes*numDim, MPI_DOUBLE, 0, 8728, mpi_info->comm, &status);
358          if ( mpi_error != MPI_SUCCESS ) {
359                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempCoords failed");
360                return NULL;
361          }
362          chunkNodes = tempInts[numNodes*3];
363    #endif
364        }   /* Worker */
365    printf("ksteube CPU=%d/%d name='%s' numDim=%d numNodes=%d chunkNodes=%d\n", mpi_info->rank, mpi_info->size, name, numDim, numNodes, chunkNodes);
366    #if 0
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    Finley_Mesh_prepare(mesh_p);      /* Copy tempMem into mesh_p */
380            Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
381        for (i0=0; i0<numNodes; i0++) {
382          mesh_p->Nodes->Id[i0]                 = tempInts[0+i0];
383          mesh_p->Nodes->globalDegreesOfFreedom[i0]     = tempInts[numNodes+i0];
384          mesh_p->Nodes->Tag[i0]                = tempInts[numNodes*2+i0];
385          for (i1=0; i1<numDim; i1++) {
386            mesh_p->Nodes->Coordinates[INDEX2(i1,i0,numDim)]    = tempCoords[i0*numDim+i1];
387          }
388            }
389    
390        TMPMEMFREE(tempInts);
391        TMPMEMFREE(tempCoords);
392    
393            return mesh_p; /* ksteube temp for debugging */
394    
395            /* read elements */
396            if (Finley_noError()) {
397      
398               fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
399               typeID=Finley_RefElement_getTypeId(element_type);
400               if (typeID==NoType) {
401                 sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s",element_type);
402                 Finley_setError(VALUE_ERROR,error_msg);
403               } else {
404                 /* read the elements */
405                 mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
406                 if (Finley_noError()) {
407                     Finley_ElementFile_allocTable(mesh_p->Elements, numEle);
408                     mesh_p->Elements->minColor=0;
409                     mesh_p->Elements->maxColor=numEle-1;
410                     if (Finley_noError()) {
411                        for (i0 = 0; i0 < numEle; i0++) {
412                          fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]);
413                          mesh_p->Elements->Color[i0]=i0;
414                          for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {
415                               fscanf(fileHandle_p, " %d",
416                                  &mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);
417                          } /* for i1 */
418                          fscanf(fileHandle_p, "\n");
419                        } /* for i0 */
420                     }
421                 }
422              }
423            }
424            /* get the face elements */
425            if (Finley_noError()) {
426                 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
427                 faceTypeID=Finley_RefElement_getTypeId(element_type);
428                 if (faceTypeID==NoType) {
429                   sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s for face elements",element_type);
430                   Finley_setError(VALUE_ERROR,error_msg);
431                 } else {
432                    mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
433                    if (Finley_noError()) {
434                       Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle);
435                       if (Finley_noError()) {
436                          mesh_p->FaceElements->minColor=0;
437                          mesh_p->FaceElements->maxColor=numEle-1;
438                          for (i0 = 0; i0 < numEle; i0++) {
439                            fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);
440                            mesh_p->FaceElements->Color[i0]=i0;
441                            for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {
442                                 fscanf(fileHandle_p, " %d",
443                                    &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);
444                            }   /* for i1 */
445                            fscanf(fileHandle_p, "\n");
446                          } /* for i0 */
447                       }
448                    }
449                 }
450            }
451            /* get the Contact face element */
452            if (Finley_noError()) {
453                 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
454                 contactTypeID=Finley_RefElement_getTypeId(element_type);
455                 if (contactTypeID==NoType) {
456                   sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type);
457                   Finley_setError(VALUE_ERROR,error_msg);
458                 } else {
459                   mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
460                   if (Finley_noError()) {
461                       Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);
462                       if (Finley_noError()) {
463                          mesh_p->ContactElements->minColor=0;
464                          mesh_p->ContactElements->maxColor=numEle-1;
465                          for (i0 = 0; i0 < numEle; i0++) {
466                            fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]);
467                            mesh_p->ContactElements->Color[i0]=i0;
468                            for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {
469                                fscanf(fileHandle_p, " %d",
470                                   &mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);
471                            }   /* for i1 */
472                            fscanf(fileHandle_p, "\n");
473                          } /* for i0 */
474                      }
475                   }
476                 }
477            }  
478            /* get the nodal element */
479            if (Finley_noError()) {
480                 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
481                 pointTypeID=Finley_RefElement_getTypeId(element_type);
482                 if (pointTypeID==NoType) {
483                   sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for points",element_type);
484                   Finley_setError(VALUE_ERROR,error_msg);
485                 }
486                 mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
487                 if (Finley_noError()) {
488                    Finley_ElementFile_allocTable(mesh_p->Points, numEle);
489                    if (Finley_noError()) {
490                       mesh_p->Points->minColor=0;
491                       mesh_p->Points->maxColor=numEle-1;
492                       for (i0 = 0; i0 < numEle; i0++) {
493                         fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);
494                         mesh_p->Points->Color[i0]=i0;
495                         for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {
496                             fscanf(fileHandle_p, " %d",
497                                &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);
498                         }  /* for i1 */
499                         fscanf(fileHandle_p, "\n");
500                       } /* for i0 */
501                    }
502                 }
503            }
504            /* get the name tags */
505            if (Finley_noError()) {
506               if (feof(fileHandle_p) == 0) {
507                  fscanf(fileHandle_p, "%s\n", name);
508                  while (feof(fileHandle_p) == 0) {
509                       fscanf(fileHandle_p, "%s %d\n", name, &tag_key);
510                       Finley_Mesh_addTagMap(mesh_p,name,tag_key);
511                  }
512               }
513            }
514         }
515         /* close file */
516         fclose(fileHandle_p);
517      
518         /*   resolve id's : */
519         /* rearrange elements: */
520      
521         if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
522         if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
523      
524         /* that's it */
525         #ifdef Finley_TRACE
526         printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);
527         #endif
528    
529    /* that's it */    /* that's it */
530    #ifdef Finley_TRACE    if (! Finley_noError()) {
531    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);  
532    }    }
533      Paso_MPIInfo_free( mpi_info );
534    return mesh_p;    return mesh_p;
535  }  }

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

  ViewVC Help
Powered by ViewVC 1.1.26