/[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 1312 by ksteube, Mon Sep 24 06:18:44 2007 UTC revision 1360 by ksteube, Thu Dec 13 05:05:00 2007 UTC
# Line 229  Finley_Mesh* Finley_Mesh_read(char* fnam Line 229  Finley_Mesh* Finley_Mesh_read(char* fnam
229    Paso_MPIInfo_free( mpi_info );    Paso_MPIInfo_free( mpi_info );
230    return mesh_p;    return mesh_p;
231  }  }
232    
233    Finley_Mesh* Finley_Mesh_read_MPI(char* fname,index_t order, index_t reduced_order,  bool_t optimize)
234    
235    {
236    
237      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      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    #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         /* 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        /* 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 */
530      if (! Finley_noError()) {
531           Finley_Mesh_free(mesh_p);
532      }
533      Paso_MPIInfo_free( mpi_info );
534      return mesh_p;
535    }

Legend:
Removed from v.1312  
changed lines
  Added in v.1360

  ViewVC Help
Powered by ViewVC 1.1.26