/[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 1725 by ksteube, Tue Aug 26 03:31:49 2008 UTC revision 1811 by ksteube, Thu Sep 25 23:11:13 2008 UTC
# Line 1  Line 1 
1    
 /* $Id$ */  
   
2  /*******************************************************  /*******************************************************
3   *  *
4   *           Copyright 2003-2007 by ACceSS MNRF  * Copyright (c) 2003-2008 by University of Queensland
5   *       Copyright 2007 by University of Queensland  * Earth Systems Science Computational Center (ESSCC)
6   *  * http://www.uq.edu.au/esscc
7   *                http://esscc.uq.edu.au  *
8   *        Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
9   *  Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
10   *     http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
11   *  *
12   *******************************************************/  *******************************************************/
13    
14    
15  /**************************************************************/  /**************************************************************/
16    
# Line 19  Line 18 
18    
19  /**************************************************************/  /**************************************************************/
20    
21    #include <ctype.h>
22  #include "Mesh.h"  #include "Mesh.h"
23    
24  /**************************************************************/  /**************************************************************/
# Line 106  Finley_Mesh* Finley_Mesh_read(char* fnam Line 106  Finley_Mesh* Finley_Mesh_read(char* fnam
106                      for (i0 = 0; i0 < numEle; i0++) {                      for (i0 = 0; i0 < numEle; i0++) {
107                        fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]);                        fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]);
108                        mesh_p->Elements->Color[i0]=i0;                        mesh_p->Elements->Color[i0]=i0;
109                          mesh_p->Elements->Owner[i0]=0;
110                        for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {                        for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {
111                             fscanf(fileHandle_p, " %d",                             fscanf(fileHandle_p, " %d",
112                                &mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);                                &mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);
# Line 133  Finley_Mesh* Finley_Mesh_read(char* fnam Line 134  Finley_Mesh* Finley_Mesh_read(char* fnam
134                        for (i0 = 0; i0 < numEle; i0++) {                        for (i0 = 0; i0 < numEle; i0++) {
135                          fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);                          fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);
136                          mesh_p->FaceElements->Color[i0]=i0;                          mesh_p->FaceElements->Color[i0]=i0;
137                            mesh_p->FaceElements->Owner[i0]=0;
138                          for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {                          for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {
139                               fscanf(fileHandle_p, " %d",                               fscanf(fileHandle_p, " %d",
140                                  &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);                                  &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);
# Line 160  Finley_Mesh* Finley_Mesh_read(char* fnam Line 162  Finley_Mesh* Finley_Mesh_read(char* fnam
162                        for (i0 = 0; i0 < numEle; i0++) {                        for (i0 = 0; i0 < numEle; i0++) {
163                          fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]);                          fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]);
164                          mesh_p->ContactElements->Color[i0]=i0;                          mesh_p->ContactElements->Color[i0]=i0;
165                            mesh_p->ContactElements->Owner[i0]=0;
166                          for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {                          for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {
167                              fscanf(fileHandle_p, " %d",                              fscanf(fileHandle_p, " %d",
168                                 &mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);                                 &mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);
# Line 187  Finley_Mesh* Finley_Mesh_read(char* fnam Line 190  Finley_Mesh* Finley_Mesh_read(char* fnam
190                     for (i0 = 0; i0 < numEle; i0++) {                     for (i0 = 0; i0 < numEle; i0++) {
191                       fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);                       fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);
192                       mesh_p->Points->Color[i0]=i0;                       mesh_p->Points->Color[i0]=i0;
193                         mesh_p->Points->Owner[i0]=0;
194                       for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {                       for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {
195                           fscanf(fileHandle_p, " %d",                           fscanf(fileHandle_p, " %d",
196                              &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);                              &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);
# Line 242  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 246  Finley_Mesh* Finley_Mesh_read_MPI(char*
246    double time0=Finley_timer();    double time0=Finley_timer();
247    FILE *fileHandle_p = NULL;    FILE *fileHandle_p = NULL;
248    ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;    ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;
249      Finley_TagMap* tag_map;
250    index_t tag_key;    index_t tag_key;
251    
252    Finley_resetError();    Finley_resetError();
# Line 308  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 313  Finley_Mesh* Finley_Mesh_read_MPI(char*
313          for (i0=0; i0<chunkSize*numDim; i0++) tempCoords[i0] = -1.0;          for (i0=0; i0<chunkSize*numDim; i0++) tempCoords[i0] = -1.0;
314          chunkNodes = 0;          chunkNodes = 0;
315          for (i1=0; i1<chunkSize; i1++) {          for (i1=0; i1<chunkSize; i1++) {
316            if (totalNodes >= numNodes) break;    /* Maybe end the infinite loop */            if (totalNodes >= numNodes) break;    /* End of inner loop */
317                if (1 == numDim)                if (1 == numDim)
318          fscanf(fileHandle_p, "%d %d %d %le\n",          fscanf(fileHandle_p, "%d %d %d %le\n",
319            &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],            &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
# Line 343  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 348  Finley_Mesh* Finley_Mesh_read_MPI(char*
348                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempCoords failed");                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempCoords failed");
349                  return NULL;                  return NULL;
350            }            }
           nextCPU++;  
351          }          }
352  #endif  #endif
353          if (totalNodes >= numNodes) break;  /* Maybe end the infinite loop */          nextCPU++;
354            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
355            if (nextCPU > mpi_info->size) break; /* End infinite loop */
356        } /* Infinite loop */        } /* Infinite loop */
357      }   /* End master */      }   /* End master */
358      else {  /* Worker */      else {  /* Worker */
# Line 368  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 374  Finley_Mesh* Finley_Mesh_read_MPI(char*
374  #endif  #endif
375      }   /* Worker */      }   /* Worker */
376    
 #if 0  
         /* Display the temp mem for debugging */  
         printf("ksteube Nodes tempInts\n");  
         for (i0=0; i0<chunkSize*3; i0++) {  
           printf(" %2d", tempInts[i0]);  
           if (i0%chunkSize==chunkSize-1) printf("\n");  
         }  
         printf("ksteube tempCoords:\n");  
         for (i0=0; i0<chunkSize*numDim; i0++) {  
           printf(" %20.15e", tempCoords[i0]);  
           if (i0%numDim==numDim-1) printf("\n");  
         }  
             printf("ksteube numDim=%d numNodes=%d mesh name='%s' chunkNodes=%d numNodes=%d\n", numDim, numNodes, name, chunkNodes, numNodes);  
 #endif  
   
377      /* Copy node data from tempMem to mesh_p */      /* Copy node data from tempMem to mesh_p */
378          Finley_NodeFile_allocTable(mesh_p->Nodes, chunkNodes);          Finley_NodeFile_allocTable(mesh_p->Nodes, chunkNodes);
379          if (Finley_noError()) {          if (Finley_noError()) {
# Line 427  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 418  Finley_Mesh* Finley_Mesh_read_MPI(char*
418            }            }
419      }      }
420    
421        /* Read the element data */        /* Allocate the ElementFile */
422        mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);        mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
423        numNodes = mesh_p->Elements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */        numNodes = mesh_p->Elements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
424    
425          /* Read the element data */
426        if (Finley_noError()) {        if (Finley_noError()) {
427      int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;      int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
428      int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */      int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
# Line 440  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 432  Finley_Mesh* Finley_Mesh_read_MPI(char*
432          for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;          for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
433          chunkEle = 0;          chunkEle = 0;
434          for (i0=0; i0<chunkSize; i0++) {          for (i0=0; i0<chunkSize; i0++) {
435            if (totalEle >= numEle) break; /* End infinite loop */            if (totalEle >= numEle) break; /* End inner loop */
436            fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);            fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
437            for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);            for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
438            fscanf(fileHandle_p, "\n");            fscanf(fileHandle_p, "\n");
# Line 457  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 449  Finley_Mesh* Finley_Mesh_read_MPI(char*
449                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Elements failed");                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Elements failed");
450                  return NULL;                  return NULL;
451            }            }
           nextCPU++;  
452          }          }
453  #endif  #endif
454          if (totalEle >= numEle) break; /* End infinite loop */          nextCPU++;
455            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
456            if (nextCPU > mpi_info->size) break; /* End infinite loop */
457        } /* Infinite loop */        } /* Infinite loop */
458      }   /* End master */      }   /* End master */
459      else {  /* Worker */      else {  /* Worker */
# Line 477  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 470  Finley_Mesh* Finley_Mesh_read_MPI(char*
470  #endif  #endif
471      }   /* Worker */      }   /* Worker */
472    
 #if 0  
     /* Display the temp mem for debugging */  
     for (i0=0; i0<chunkSize*(numNodes+2); i0++) {  
       printf(" %2d", tempInts[i0]);  
       if (i0%(numNodes+2)==numNodes+2-1) printf("\n");  
     }  
 #endif  
   
473      /* Copy Element data from tempInts to mesh_p */      /* Copy Element data from tempInts to mesh_p */
474      Finley_ElementFile_allocTable(mesh_p->Elements, chunkEle);      Finley_ElementFile_allocTable(mesh_p->Elements, chunkEle);
475          mesh_p->Elements->minColor=0;          mesh_p->Elements->minColor=0;
# Line 494  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 479  Finley_Mesh* Finley_Mesh_read_MPI(char*
479        for (i0=0; i0<chunkEle; i0++) {        for (i0=0; i0<chunkEle; i0++) {
480          mesh_p->Elements->Id[i0]    = tempInts[i0*(2+numNodes)+0];          mesh_p->Elements->Id[i0]    = tempInts[i0*(2+numNodes)+0];
481          mesh_p->Elements->Tag[i0]   = tempInts[i0*(2+numNodes)+1];          mesh_p->Elements->Tag[i0]   = tempInts[i0*(2+numNodes)+1];
482                mesh_p->Elements->Owner[i0]  =mpi_info->rank;
483              mesh_p->Elements->Color[i0] = i0;              mesh_p->Elements->Color[i0] = i0;
484          for (i1 = 0; i1 < numNodes; i1++) {          for (i1 = 0; i1 < numNodes; i1++) {
485            mesh_p->Elements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];            mesh_p->Elements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
# Line 502  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 488  Finley_Mesh* Finley_Mesh_read_MPI(char*
488          }          }
489    
490      TMPMEMFREE(tempInts);      TMPMEMFREE(tempInts);
491          } /* end of Read the element data */
492    
493            /* read face elements */
494    
495        /* Read the element typeID */
496            if (Finley_noError()) {
497          if (mpi_info->rank == 0) {
498                fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
499                typeID=Finley_RefElement_getTypeId(element_type);
500          }
501    #ifdef PASO_MPI
502          if (mpi_info->size > 1) {
503            int temp1[2], mpi_error;
504            temp1[0] = (int) typeID;
505            temp1[1] = numEle;
506            mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
507            if (mpi_error != MPI_SUCCESS) {
508              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
509              return NULL;
510            }
511            typeID = (ElementTypeId) temp1[0];
512            numEle = temp1[1];
513          }
514    #endif
515              if (typeID==NoType) {
516                sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
517                Finley_setError(VALUE_ERROR, error_msg);
518              }
519        }
520    
521          /* Allocate the ElementFile */
522          mesh_p->FaceElements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
523          numNodes = mesh_p->FaceElements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
524    
525          /* Read the face element data */
526          if (Finley_noError()) {
527        int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
528        int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
529        /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
530        if (mpi_info->rank == 0) {  /* Master */
531          for (;;) {            /* Infinite loop */
532            for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
533            chunkEle = 0;
534            for (i0=0; i0<chunkSize; i0++) {
535              if (totalEle >= numEle) break; /* End inner loop */
536              fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
537              for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
538              fscanf(fileHandle_p, "\n");
539              totalEle++;
540              chunkEle++;
541            }
542    #ifdef PASO_MPI
543            /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
544            if (nextCPU < mpi_info->size) {
545                  int mpi_error;
546              tempInts[chunkSize*(2+numNodes)] = chunkEle;
547              mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81723, mpi_info->comm);
548              if ( mpi_error != MPI_SUCCESS ) {
549                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for FaceElements failed");
550                    return NULL;
551              }
552            }
553    #endif
554            nextCPU++;
555            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
556            if (nextCPU > mpi_info->size) break; /* End infinite loop */
557          } /* Infinite loop */
558        }   /* End master */
559        else {  /* Worker */
560    #ifdef PASO_MPI
561          /* Each worker receives one message */
562          MPI_Status status;
563          int mpi_error;
564          mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81723, mpi_info->comm, &status);
565          if ( mpi_error != MPI_SUCCESS ) {
566                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for FaceElements failed");
567                return NULL;
568          }
569          chunkEle = tempInts[chunkSize*(2+numNodes)];
570    #endif
571        }   /* Worker */
572    
573        /* Copy Element data from tempInts to mesh_p */
574        Finley_ElementFile_allocTable(mesh_p->FaceElements, chunkEle);
575            mesh_p->FaceElements->minColor=0;
576            mesh_p->FaceElements->maxColor=chunkEle-1;
577            if (Finley_noError()) {
578              #pragma omp parallel for private (i0, i1)
579          for (i0=0; i0<chunkEle; i0++) {
580            mesh_p->FaceElements->Id[i0]    = tempInts[i0*(2+numNodes)+0];
581            mesh_p->FaceElements->Tag[i0]   = tempInts[i0*(2+numNodes)+1];
582                mesh_p->FaceElements->Owner[i0]  =mpi_info->rank;
583                mesh_p->FaceElements->Color[i0] = i0;
584            for (i1 = 0; i1 < numNodes; i1++) {
585              mesh_p->FaceElements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
586            }
587              }
588            }
589    
590        TMPMEMFREE(tempInts);
591          } /* end of Read the face element data */
592    
593            /* read contact elements */
594    
595        /* Read the element typeID */
596            if (Finley_noError()) {
597          if (mpi_info->rank == 0) {
598                fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
599                typeID=Finley_RefElement_getTypeId(element_type);
600          }
601    #ifdef PASO_MPI
602          if (mpi_info->size > 1) {
603            int temp1[2], mpi_error;
604            temp1[0] = (int) typeID;
605            temp1[1] = numEle;
606            mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
607            if (mpi_error != MPI_SUCCESS) {
608              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
609              return NULL;
610            }
611            typeID = (ElementTypeId) temp1[0];
612            numEle = temp1[1];
613          }
614    #endif
615              if (typeID==NoType) {
616                sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
617                Finley_setError(VALUE_ERROR, error_msg);
618              }
619        }
620    
621          /* Allocate the ElementFile */
622          mesh_p->ContactElements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
623          numNodes = mesh_p->ContactElements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
624    
625          /* Read the contact element data */
626          if (Finley_noError()) {
627        int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
628        int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
629        /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
630        if (mpi_info->rank == 0) {  /* Master */
631          for (;;) {            /* Infinite loop */
632            for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
633            chunkEle = 0;
634            for (i0=0; i0<chunkSize; i0++) {
635              if (totalEle >= numEle) break; /* End inner loop */
636              fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
637              for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
638              fscanf(fileHandle_p, "\n");
639              totalEle++;
640              chunkEle++;
641            }
642    #ifdef PASO_MPI
643            /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
644            if (nextCPU < mpi_info->size) {
645                  int mpi_error;
646              tempInts[chunkSize*(2+numNodes)] = chunkEle;
647              mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81724, mpi_info->comm);
648              if ( mpi_error != MPI_SUCCESS ) {
649                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for ContactElements failed");
650                    return NULL;
651              }
652            }
653    #endif
654            nextCPU++;
655            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
656            if (nextCPU > mpi_info->size) break; /* End infinite loop */
657          } /* Infinite loop */
658        }   /* End master */
659        else {  /* Worker */
660    #ifdef PASO_MPI
661          /* Each worker receives one message */
662          MPI_Status status;
663          int mpi_error;
664          mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81724, mpi_info->comm, &status);
665          if ( mpi_error != MPI_SUCCESS ) {
666                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for ContactElements failed");
667                return NULL;
668          }
669          chunkEle = tempInts[chunkSize*(2+numNodes)];
670    #endif
671        }   /* Worker */
672    
673        /* Copy Element data from tempInts to mesh_p */
674        Finley_ElementFile_allocTable(mesh_p->ContactElements, chunkEle);
675            mesh_p->ContactElements->minColor=0;
676            mesh_p->ContactElements->maxColor=chunkEle-1;
677            if (Finley_noError()) {
678              #pragma omp parallel for private (i0, i1)
679          for (i0=0; i0<chunkEle; i0++) {
680            mesh_p->ContactElements->Id[i0] = tempInts[i0*(2+numNodes)+0];
681            mesh_p->ContactElements->Tag[i0]    = tempInts[i0*(2+numNodes)+1];
682                mesh_p->ContactElements->Owner[i0]  =mpi_info->rank;
683                mesh_p->ContactElements->Color[i0] = i0;
684            for (i1 = 0; i1 < numNodes; i1++) {
685              mesh_p->ContactElements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
686            }
687              }
688            }
689    
690        TMPMEMFREE(tempInts);
691          } /* end of Read the contact element data */
692    
693            /* read nodal elements */
694    
695        /* Read the element typeID */
696            if (Finley_noError()) {
697          if (mpi_info->rank == 0) {
698                fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
699                typeID=Finley_RefElement_getTypeId(element_type);
700          }
701    #ifdef PASO_MPI
702          if (mpi_info->size > 1) {
703            int temp1[2], mpi_error;
704            temp1[0] = (int) typeID;
705            temp1[1] = numEle;
706            mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
707            if (mpi_error != MPI_SUCCESS) {
708              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
709              return NULL;
710            }
711            typeID = (ElementTypeId) temp1[0];
712            numEle = temp1[1];
713          }
714    #endif
715              if (typeID==NoType) {
716                sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
717                Finley_setError(VALUE_ERROR, error_msg);
718              }
719        }
720    
721          /* Allocate the ElementFile */
722          mesh_p->Points=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
723          numNodes = mesh_p->Points->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
724    
725          /* Read the nodal element data */
726          if (Finley_noError()) {
727        int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
728        int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
729        /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
730        if (mpi_info->rank == 0) {  /* Master */
731          for (;;) {            /* Infinite loop */
732            for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
733            chunkEle = 0;
734            for (i0=0; i0<chunkSize; i0++) {
735              if (totalEle >= numEle) break; /* End inner loop */
736              fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
737              for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
738              fscanf(fileHandle_p, "\n");
739              totalEle++;
740              chunkEle++;
741            }
742    #ifdef PASO_MPI
743            /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
744            if (nextCPU < mpi_info->size) {
745                  int mpi_error;
746              tempInts[chunkSize*(2+numNodes)] = chunkEle;
747              mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81725, mpi_info->comm);
748              if ( mpi_error != MPI_SUCCESS ) {
749                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Points failed");
750                    return NULL;
751              }
752            }
753    #endif
754            nextCPU++;
755            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
756            if (nextCPU > mpi_info->size) break; /* End infinite loop */
757          } /* Infinite loop */
758        }   /* End master */
759        else {  /* Worker */
760    #ifdef PASO_MPI
761          /* Each worker receives one message */
762          MPI_Status status;
763          int mpi_error;
764          mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81725, mpi_info->comm, &status);
765          if ( mpi_error != MPI_SUCCESS ) {
766                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for Points failed");
767                return NULL;
768          }
769          chunkEle = tempInts[chunkSize*(2+numNodes)];
770    #endif
771        }   /* Worker */
772    
773        /* Copy Element data from tempInts to mesh_p */
774        Finley_ElementFile_allocTable(mesh_p->Points, chunkEle);
775            mesh_p->Points->minColor=0;
776            mesh_p->Points->maxColor=chunkEle-1;
777            if (Finley_noError()) {
778              #pragma omp parallel for private (i0, i1)
779          for (i0=0; i0<chunkEle; i0++) {
780            mesh_p->Points->Id[i0]  = tempInts[i0*(2+numNodes)+0];
781            mesh_p->Points->Tag[i0] = tempInts[i0*(2+numNodes)+1];
782                mesh_p->Points->Owner[i0]  =mpi_info->rank;
783                mesh_p->Points->Color[i0] = i0;
784            for (i1 = 0; i1 < numNodes; i1++) {
785              mesh_p->Points->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
786            }
787              }
788            }
789    
790        TMPMEMFREE(tempInts);
791          } /* end of Read the nodal element data */
792    
793          /* get the name tags */
794          if (Finley_noError()) {
795            char *remainder, *ptr;
796            int tag_key, len, error_code;
797            long cur_pos, end_pos;
798            if (mpi_info->rank == 0) {  /* Master */
799          /* Read the word 'Tag' */
800          if (! feof(fileHandle_p)) fscanf(fileHandle_p, "%s\n", name);
801          /* Read rest of file in one chunk, after using seek to find length */
802              cur_pos = ftell(fileHandle_p);
803              fseek(fileHandle_p, 0L, SEEK_END);
804              end_pos = ftell(fileHandle_p);
805              fseek(fileHandle_p, (long)cur_pos, SEEK_SET);
806          remainder = TMPMEMALLOC(end_pos-cur_pos+1, char);
807          if (! feof(fileHandle_p)) fread(remainder, (size_t) end_pos-cur_pos, sizeof(char), fileHandle_p);
808          remainder[end_pos-cur_pos] = 0;
809          len = strlen(remainder);    
810          while ((--len)>0 && isspace(remainder[len])) remainder[len]=0;
811          len = strlen(remainder);
812            }
813    #ifdef PASO_MPI
814            error_code = MPI_Bcast (&len, 1, MPI_INT,  0, mpi_info->comm);
815            if (error_code != MPI_SUCCESS) {
816              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of tag len failed");
817              return NULL;
818            }
819        if (mpi_info->rank != 0) {
820          remainder = TMPMEMALLOC(len+1, char);
821          remainder[0] = 0;
822        }
823            error_code = MPI_Bcast (remainder, len+1, MPI_CHAR,  0, mpi_info->comm);
824            if (error_code != MPI_SUCCESS) {
825              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of tags failed");
826              return NULL;
827            }
828    #endif
829        if (remainder[0]) {
830              ptr = remainder;
831              do {
832                sscanf(ptr, "%s %d\n", name, &tag_key);
833                if (*name) Finley_Mesh_addTagMap(mesh_p,name,tag_key);
834                ptr++;
835              } while(NULL != (ptr = strchr(ptr, '\n')) && *ptr);
836              TMPMEMFREE(remainder);
837        }
838        }        }
839    
840       }       }
841    
842       /* close file */       /* close file */
# Line 511  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 845  Finley_Mesh* Finley_Mesh_read_MPI(char*
845       /*   resolve id's : */       /*   resolve id's : */
846       /* rearrange elements: */       /* rearrange elements: */
847    
      return mesh_p; /* ksteube temp return for debugging */  
   
848       if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);       if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
849       if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);       if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
850    

Legend:
Removed from v.1725  
changed lines
  Added in v.1811

  ViewVC Help
Powered by ViewVC 1.1.26