/[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 1628 by phornby, Fri Jul 11 13:12:46 2008 UTC revision 1732 by ksteube, Wed Aug 27 05:47:03 2008 UTC
# Line 241  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 241  Finley_Mesh* Finley_Mesh_read_MPI(char*
241    char error_msg[LenErrorMsg_MAX];    char error_msg[LenErrorMsg_MAX];
242    double time0=Finley_timer();    double time0=Finley_timer();
243    FILE *fileHandle_p = NULL;    FILE *fileHandle_p = NULL;
244    ElementTypeId typeID;    ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;
   
   /* these are in unimplemented code below */  
 #if 0  
245    index_t tag_key;    index_t tag_key;
   ElementTypeId faceTypeID, contactTypeID, pointTypeID;  
 #endif  
246    
247    Finley_resetError();    Finley_resetError();
248    
# Line 270  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 265  Finley_Mesh* Finley_Mesh_read_MPI(char*
265    }    }
266    
267  #ifdef PASO_MPI  #ifdef PASO_MPI
268    /* MPI Broadcast numDim, numNodes, name */    /* MPI Broadcast numDim, numNodes, name if there are multiple MPI procs*/
269    if (mpi_info->size > 0) {    if (mpi_info->size > 1) {
270      int temp1[3], error_code;      int temp1[3], error_code;
271      temp1[0] = numDim;      temp1[0] = numDim;
272      temp1[1] = numNodes;      temp1[1] = numNodes;
# Line 294  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 289  Finley_Mesh* Finley_Mesh_read_MPI(char*
289       /* allocate mesh */       /* allocate mesh */
290       mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);       mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
291       if (Finley_noError()) {       if (Finley_noError()) {
292        /* Each CPU will get at most chunkSize nodes so the message has to be sufficiently large */
293      int chunkSize = numNodes / mpi_info->size + 1, totalNodes=0, chunkNodes=0, chunkEle=0, nextCPU=1;      int chunkSize = numNodes / mpi_info->size + 1, totalNodes=0, chunkNodes=0, chunkEle=0, nextCPU=1;
294  #ifdef PASO_MPI      int *tempInts = TMPMEMALLOC(chunkSize*3+1, index_t);        /* Stores the integer message data */
295      int mpi_error;      double *tempCoords = TMPMEMALLOC(chunkSize*numDim, double); /* Stores the double message data */
 #endif  
     int *tempInts = TMPMEMALLOC(numNodes*3+1, index_t);  
     double *tempCoords = TMPMEMALLOC(numNodes*numDim, double);  
296    
297      /*      /*
298        Read a chunk of nodes, send to worker CPU if available, copy chunk into local mesh_p        Read chunkSize nodes, send it in a chunk to worker CPU which copies chunk into its local mesh_p
299        It doesn't matter that a CPU has the wrong nodes for its elements, this is sorted out later        It doesn't matter that a CPU has the wrong nodes for its elements, this is sorted out later
300        First chunk sent to CPU 1, second to CPU 2, ...        First chunk sent to CPU 1, second to CPU 2, ...
301        Last chunk stays on CPU 0 (the master)        Last chunk stays on CPU 0 (the master)
# Line 311  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 304  Finley_Mesh* Finley_Mesh_read_MPI(char*
304    
305      if (mpi_info->rank == 0) {  /* Master */      if (mpi_info->rank == 0) {  /* Master */
306        for (;;) {            /* Infinite loop */        for (;;) {            /* Infinite loop */
307            for (i0=0; i0<chunkSize*3+1; i0++) tempInts[i0] = -1;
308            for (i0=0; i0<chunkSize*numDim; i0++) tempCoords[i0] = -1.0;
309          chunkNodes = 0;          chunkNodes = 0;
         for (i0=0; i0<numNodes*3+1; i0++) tempInts[i0] = -1;  
         for (i0=0; i0<numNodes*numDim; i0++) tempCoords[i0] = -1.0;  
310          for (i1=0; i1<chunkSize; i1++) {          for (i1=0; i1<chunkSize; i1++) {
311            if (totalNodes >= numNodes) break;            if (totalNodes >= numNodes) break;    /* End of inner loop */
312                if (1 == numDim)                if (1 == numDim)
313          fscanf(fileHandle_p, "%d %d %d %le\n",          fscanf(fileHandle_p, "%d %d %d %le\n",
314            &tempInts[0+i1], &tempInts[numNodes+i1], &tempInts[numNodes*2+i1],            &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
315            &tempCoords[i1*numDim+0]);            &tempCoords[i1*numDim+0]);
316                if (2 == numDim)                if (2 == numDim)
317          fscanf(fileHandle_p, "%d %d %d %le %le\n",          fscanf(fileHandle_p, "%d %d %d %le %le\n",
318            &tempInts[0+i1], &tempInts[numNodes+i1], &tempInts[numNodes*2+i1],            &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
319            &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1]);            &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1]);
320                if (3 == numDim)                if (3 == numDim)
321          fscanf(fileHandle_p, "%d %d %d %le %le %le\n",          fscanf(fileHandle_p, "%d %d %d %le %le %le\n",
322            &tempInts[0+i1], &tempInts[numNodes+i1], &tempInts[numNodes*2+i1],            &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
323            &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1], &tempCoords[i1*numDim+2]);            &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1], &tempCoords[i1*numDim+2]);
324            totalNodes++;            totalNodes++; /* When do we quit the infinite loop? */
325            chunkNodes++;            chunkNodes++; /* How many nodes do we actually have in this chunk? It may be smaller than chunkSize. */
326            }
327            if (chunkNodes > chunkSize) {
328                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: error reading chunks of mesh, data too large for message size");
329                  return NULL;
330          }          }
         /* Eventually we'll send chunk of nodes to each CPU numbered 1 ... mpi_info->size-1, here goes one of them */  
         if (nextCPU < mpi_info->size) {  
331  #ifdef PASO_MPI  #ifdef PASO_MPI
332            tempInts[numNodes*3] = chunkNodes;          /* Eventually we'll send chunkSize nodes to each CPU numbered 1 ... mpi_info->size-1, here goes one of them */
333            /* ksteube The size of this message can and should be brought down to chunkNodes*3+1, must re-org tempInts */          if (nextCPU < mpi_info->size) {
334            mpi_error = MPI_Send(tempInts, numNodes*3+1, MPI_INT, nextCPU, 81720, mpi_info->comm);                int mpi_error;
335              tempInts[chunkSize*3] = chunkNodes;   /* The message has one more int to send chunkNodes */
336              mpi_error = MPI_Send(tempInts, chunkSize*3+1, MPI_INT, nextCPU, 81720, mpi_info->comm);
337            if ( mpi_error != MPI_SUCCESS ) {            if ( mpi_error != MPI_SUCCESS ) {
338                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts failed");                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts failed");
339                  return NULL;                  return NULL;
340            }            }
341            mpi_error = MPI_Send(tempCoords, numNodes*numDim, MPI_DOUBLE, nextCPU, 81721, mpi_info->comm);            mpi_error = MPI_Send(tempCoords, chunkSize*numDim, MPI_DOUBLE, nextCPU, 81721, mpi_info->comm);
342            if ( mpi_error != MPI_SUCCESS ) {            if ( mpi_error != MPI_SUCCESS ) {
343                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempCoords failed");                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempCoords failed");
344                  return NULL;                  return NULL;
345            }            }
 #endif  
           nextCPU++;  
346          }          }
347          if (totalNodes >= numNodes) break;  #endif
348            nextCPU++;
349            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
350            if (nextCPU > mpi_info->size) break; /* End infinite loop */
351        } /* Infinite loop */        } /* Infinite loop */
352      }   /* End master */      }   /* End master */
353      else {  /* Worker */      else {  /* Worker */
354  #ifdef PASO_MPI  #ifdef PASO_MPI
355        /* Each worker receives two messages */        /* Each worker receives two messages */
356        MPI_Status status;        MPI_Status status;
357        mpi_error = MPI_Recv(tempInts, numNodes*3+1, MPI_INT, 0, 81720, mpi_info->comm, &status);            int mpi_error;
358          mpi_error = MPI_Recv(tempInts, chunkSize*3+1, MPI_INT, 0, 81720, mpi_info->comm, &status);
359        if ( mpi_error != MPI_SUCCESS ) {        if ( mpi_error != MPI_SUCCESS ) {
360              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts failed");              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts failed");
361              return NULL;              return NULL;
362        }        }
363        mpi_error = MPI_Recv(tempCoords, numNodes*numDim, MPI_DOUBLE, 0, 81721, mpi_info->comm, &status);        mpi_error = MPI_Recv(tempCoords, chunkSize*numDim, MPI_DOUBLE, 0, 81721, mpi_info->comm, &status);
364        if ( mpi_error != MPI_SUCCESS ) {        if ( mpi_error != MPI_SUCCESS ) {
365              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempCoords failed");              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempCoords failed");
366              return NULL;              return NULL;
367        }        }
368        chunkNodes = tempInts[numNodes*3];        chunkNodes = tempInts[chunkSize*3];   /* How many nodes are in this workers chunk? */
369  #endif  #endif
370      }   /* Worker */      }   /* Worker */
371    
 #if 0  
         /* Display the temp mem for debugging */  
         printf("ksteube tempInts totalNodes=%d:\n", totalNodes);  
         for (i0=0; i0<numNodes*3; i0++) {  
           printf(" %2d", tempInts[i0]);  
           if (i0%numNodes==numNodes-1) printf("\n");  
         }  
         printf("ksteube tempCoords:\n");  
         for (i0=0; i0<chunkNodes*numDim; i0++) {  
           printf(" %20.15e", tempCoords[i0]);  
           if (i0%numDim==numDim-1) printf("\n");  
         }  
 #endif  
   
     printf("ksteube chunkNodes=%d numNodes=%d\n", chunkNodes, numNodes);  
372      /* Copy node data from tempMem to mesh_p */      /* Copy node data from tempMem to mesh_p */
373          Finley_NodeFile_allocTable(mesh_p->Nodes, chunkNodes);          Finley_NodeFile_allocTable(mesh_p->Nodes, chunkNodes);
374          if (Finley_noError()) {          if (Finley_noError()) {
375        for (i0=0; i0<chunkNodes; i0++) {        for (i0=0; i0<chunkNodes; i0++) {
376          mesh_p->Nodes->Id[i0]               = tempInts[0+i0];          mesh_p->Nodes->Id[i0]               = tempInts[0+i0];
377          mesh_p->Nodes->globalDegreesOfFreedom[i0]       = tempInts[numNodes+i0];          mesh_p->Nodes->globalDegreesOfFreedom[i0]       = tempInts[chunkSize+i0];
378          mesh_p->Nodes->Tag[i0]              = tempInts[numNodes*2+i0];          mesh_p->Nodes->Tag[i0]              = tempInts[chunkSize*2+i0];
379          for (i1=0; i1<numDim; i1++) {          for (i1=0; i1<numDim; i1++) {
380            mesh_p->Nodes->Coordinates[INDEX2(i1,i0,numDim)]  = tempCoords[i0*numDim+i1];            mesh_p->Nodes->Coordinates[INDEX2(i1,i0,numDim)]  = tempCoords[i0*numDim+i1];
381          }          }
# Line 410  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 394  Finley_Mesh* Finley_Mesh_read_MPI(char*
394              typeID=Finley_RefElement_getTypeId(element_type);              typeID=Finley_RefElement_getTypeId(element_type);
395        }        }
396  #ifdef PASO_MPI  #ifdef PASO_MPI
397        if (mpi_info->size > 0) {        if (mpi_info->size > 1) {
398          int temp1[3];          int temp1[2], mpi_error;
399          temp1[0] = typeID;          temp1[0] = (int) typeID;
400          temp1[1] = numEle;          temp1[1] = numEle;
401          mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);          mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
402          if (mpi_error != MPI_SUCCESS) {          if (mpi_error != MPI_SUCCESS) {
403            Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");            Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
404            return NULL;            return NULL;
405          }          }
406          typeID = temp1[0];          typeID = (ElementTypeId) temp1[0];
407          numEle = temp1[1];          numEle = temp1[1];
408        }        }
409  #endif  #endif
# Line 429  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 413  Finley_Mesh* Finley_Mesh_read_MPI(char*
413            }            }
414      }      }
415    
416        /* Read the element data */        /* Allocate the ElementFile */
417        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);
418        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 */
419    
420          /* Read the element data */
421        if (Finley_noError()) {        if (Finley_noError()) {
422      int *tempInts = TMPMEMALLOC(numEle*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */      int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
423      int chunkSize = numEle / mpi_info->size, totalEle=0, nextCPU=1;      int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
424  #ifdef PASO_MPI      /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
     int mpi_error;  
 #endif  
     if (numEle % mpi_info->size != 0) chunkSize++; /* Remainder from numEle / mpi_info->size will be spread out one-per-CPU */  
425      if (mpi_info->rank == 0) {  /* Master */      if (mpi_info->rank == 0) {  /* Master */
426        for (;;) {            /* Infinite loop */        for (;;) {            /* Infinite loop */
427            for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
428          chunkEle = 0;          chunkEle = 0;
         for (i0=0; i0<numEle*(2+numNodes)+1; i0++) tempInts[i0] = -1;  
429          for (i0=0; i0<chunkSize; i0++) {          for (i0=0; i0<chunkSize; i0++) {
430            if (totalEle >= numEle) break; /* End infinite loop */            if (totalEle >= numEle) break; /* End inner loop */
431            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]);
432            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]);
433            fscanf(fileHandle_p, "\n");            fscanf(fileHandle_p, "\n");
434            totalEle++;            totalEle++;
435            chunkEle++;            chunkEle++;
436          }          }
         /* Eventually we'll send chunk of nodes to each CPU except 0 itself, here goes one of them */  
         if (nextCPU < mpi_info->size) {  
437  #ifdef PASO_MPI  #ifdef PASO_MPI
438            tempInts[numEle*(2+numNodes)] = chunkEle;          /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
439  printf("ksteube CPU=%d/%d send to %d\n", mpi_info->rank, mpi_info->size, nextCPU);          if (nextCPU < mpi_info->size) {
440            mpi_error = MPI_Send(tempInts, numEle*(2+numNodes)+1, MPI_INT, nextCPU, 81722, mpi_info->comm);                int mpi_error;
441              tempInts[chunkSize*(2+numNodes)] = chunkEle;
442              mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81722, mpi_info->comm);
443            if ( mpi_error != MPI_SUCCESS ) {            if ( mpi_error != MPI_SUCCESS ) {
444                  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");
445                  return NULL;                  return NULL;
446            }            }
 #endif  
           nextCPU++;  
447          }          }
448          if (totalEle >= numEle) break; /* End infinite loop */  #endif
449            nextCPU++;
450            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
451            if (nextCPU > mpi_info->size) break; /* End infinite loop */
452        } /* Infinite loop */        } /* Infinite loop */
453      }   /* End master */      }   /* End master */
454      else {  /* Worker */      else {  /* Worker */
455  #ifdef PASO_MPI  #ifdef PASO_MPI
456        /* Each worker receives two messages */        /* Each worker receives one message */
457        MPI_Status status;        MPI_Status status;
458  printf("ksteube CPU=%d/%d recv on %d\n", mpi_info->rank, mpi_info->size, mpi_info->rank);        int mpi_error;
459        mpi_error = MPI_Recv(tempInts, numEle*(2+numNodes)+1, MPI_INT, 0, 81722, mpi_info->comm, &status);        mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81722, mpi_info->comm, &status);
460        if ( mpi_error != MPI_SUCCESS ) {        if ( mpi_error != MPI_SUCCESS ) {
461              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for Elements failed");              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for Elements failed");
462              return NULL;              return NULL;
463        }        }
464        chunkEle = tempInts[numEle*(2+numNodes)];        chunkEle = tempInts[chunkSize*(2+numNodes)];
465  #endif  #endif
466      }   /* Worker */      }   /* Worker */
 #if 1  
     /* Display the temp mem for debugging */  
     printf("ksteube tempInts numEle=%d chunkEle=%d AAA:\n", numEle, chunkEle);  
     for (i0=0; i0<numEle*(numNodes+2); i0++) {  
       printf(" %2d", tempInts[i0]);  
       if (i0%(numNodes+2)==numNodes+2-1) printf("\n");  
     }  
 #endif  
467    
468      /* Copy Element data from tempInts to mesh_p */      /* Copy Element data from tempInts to mesh_p */
469      Finley_ElementFile_allocTable(mesh_p->Elements, chunkEle);      Finley_ElementFile_allocTable(mesh_p->Elements, chunkEle);
# Line 499  printf("ksteube CPU=%d/%d recv on %d\n", Line 474  printf("ksteube CPU=%d/%d recv on %d\n",
474        for (i0=0; i0<chunkEle; i0++) {        for (i0=0; i0<chunkEle; i0++) {
475          mesh_p->Elements->Id[i0]    = tempInts[i0*(2+numNodes)+0];          mesh_p->Elements->Id[i0]    = tempInts[i0*(2+numNodes)+0];
476          mesh_p->Elements->Tag[i0]   = tempInts[i0*(2+numNodes)+1];          mesh_p->Elements->Tag[i0]   = tempInts[i0*(2+numNodes)+1];
477                mesh_p->Elements->Color[i0] = i0;
478          for (i1 = 0; i1 < numNodes; i1++) {          for (i1 = 0; i1 < numNodes; i1++) {
479            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];
480          }          }
# Line 506  printf("ksteube CPU=%d/%d recv on %d\n", Line 482  printf("ksteube CPU=%d/%d recv on %d\n",
482          }          }
483    
484      TMPMEMFREE(tempInts);      TMPMEMFREE(tempInts);
485        }        } /* end of Read the element data */
486    
487            /* read face elements */
488    
489  #if 0 /* this is the original code for reading elements */      /* Read the element typeID */
         /* read elements */  
490          if (Finley_noError()) {          if (Finley_noError()) {
491          if (mpi_info->rank == 0) {
492             fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);              fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
493             typeID=Finley_RefElement_getTypeId(element_type);              typeID=Finley_RefElement_getTypeId(element_type);
494             if (typeID==NoType) {        }
495               sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s",element_type);  #ifdef PASO_MPI
496               Finley_setError(VALUE_ERROR,error_msg);        if (mpi_info->size > 1) {
497             } else {          int temp1[2], mpi_error;
498               /* read the elements */          temp1[0] = (int) typeID;
499               mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);          temp1[1] = numEle;
500               if (Finley_noError()) {          mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
501                   Finley_ElementFile_allocTable(mesh_p->Elements, numEle);          if (mpi_error != MPI_SUCCESS) {
502                   mesh_p->Elements->minColor=0;            Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
503                   mesh_p->Elements->maxColor=numEle-1;            return NULL;
504                   if (Finley_noError()) {          }
505                      for (i0 = 0; i0 < numEle; i0++) {          typeID = (ElementTypeId) temp1[0];
506                        fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]);          numEle = temp1[1];
507                        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 */  
                  }  
              }  
           }  
         }  
508  #endif  #endif
509              if (typeID==NoType) {
510                sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
511                Finley_setError(VALUE_ERROR, error_msg);
512              }
513        }
514    
515  printf("ksteube CPU=%d/%d Element typeID=%d\n", mpi_info->rank, mpi_info->size, typeID);        /* Allocate the ElementFile */
516          mesh_p->FaceElements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
517          numNodes = mesh_p->FaceElements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
518    
519  #if 1        /* Read the face element data */
520          if (Finley_noError()) {
521        int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
522        int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
523        /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
524        if (mpi_info->rank == 0) {  /* Master */
525          for (;;) {            /* Infinite loop */
526            for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
527            chunkEle = 0;
528            for (i0=0; i0<chunkSize; i0++) {
529              if (totalEle >= numEle) break; /* End inner loop */
530              fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
531              for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
532              fscanf(fileHandle_p, "\n");
533              totalEle++;
534              chunkEle++;
535            }
536    #ifdef PASO_MPI
537            /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
538            if (nextCPU < mpi_info->size) {
539                  int mpi_error;
540              tempInts[chunkSize*(2+numNodes)] = chunkEle;
541              mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81723, mpi_info->comm);
542              if ( mpi_error != MPI_SUCCESS ) {
543                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for FaceElements failed");
544                    return NULL;
545              }
546            }
547    #endif
548            nextCPU++;
549            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
550            if (nextCPU > mpi_info->size) break; /* End infinite loop */
551          } /* Infinite loop */
552        }   /* End master */
553        else {  /* Worker */
554    #ifdef PASO_MPI
555          /* Each worker receives one message */
556          MPI_Status status;
557          int mpi_error;
558          mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81723, mpi_info->comm, &status);
559          if ( mpi_error != MPI_SUCCESS ) {
560                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for FaceElements failed");
561                return NULL;
562          }
563          chunkEle = tempInts[chunkSize*(2+numNodes)];
564    #endif
565        }   /* Worker */
566    
567    /* Define other structures to keep mesh_write from crashing */      /* Copy Element data from tempInts to mesh_p */
568        Finley_ElementFile_allocTable(mesh_p->FaceElements, chunkEle);
569            mesh_p->FaceElements->minColor=0;
570            mesh_p->FaceElements->maxColor=chunkEle-1;
571            if (Finley_noError()) {
572              #pragma omp parallel for private (i0, i1)
573          for (i0=0; i0<chunkEle; i0++) {
574            mesh_p->FaceElements->Id[i0]    = tempInts[i0*(2+numNodes)+0];
575            mesh_p->FaceElements->Tag[i0]   = tempInts[i0*(2+numNodes)+1];
576                mesh_p->FaceElements->Color[i0] = i0;
577            for (i1 = 0; i1 < numNodes; i1++) {
578              mesh_p->FaceElements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
579            }
580              }
581            }
582    
583    mesh_p->FaceElements=Finley_ElementFile_alloc(0, mesh_p->order, mesh_p->reduced_order, mpi_info);      TMPMEMFREE(tempInts);
584    Finley_ElementFile_allocTable(mesh_p->FaceElements, 0);        } /* end of Read the face element data */
585    
586    mesh_p->ContactElements=Finley_ElementFile_alloc(0, mesh_p->order, mesh_p->reduced_order, mpi_info);          /* read contact elements */
   Finley_ElementFile_allocTable(mesh_p->ContactElements, 0);  
587    
588    mesh_p->Points=Finley_ElementFile_alloc(0, mesh_p->order, mesh_p->reduced_order, mpi_info);      /* Read the element typeID */
589    Finley_ElementFile_allocTable(mesh_p->Points, 0);          if (Finley_noError()) {
590          if (mpi_info->rank == 0) {
591                fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
592                typeID=Finley_RefElement_getTypeId(element_type);
593          }
594    #ifdef PASO_MPI
595          if (mpi_info->size > 1) {
596            int temp1[2], mpi_error;
597            temp1[0] = (int) typeID;
598            temp1[1] = numEle;
599            mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
600            if (mpi_error != MPI_SUCCESS) {
601              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
602              return NULL;
603            }
604            typeID = (ElementTypeId) temp1[0];
605            numEle = temp1[1];
606          }
607    #endif
608              if (typeID==NoType) {
609                sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
610                Finley_setError(VALUE_ERROR, error_msg);
611              }
612        }
613    
614          /* Allocate the ElementFile */
615          mesh_p->ContactElements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
616          numNodes = mesh_p->ContactElements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
617    
618          /* Read the contact element data */
619          if (Finley_noError()) {
620        int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
621        int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
622        /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
623        if (mpi_info->rank == 0) {  /* Master */
624          for (;;) {            /* Infinite loop */
625            for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
626            chunkEle = 0;
627            for (i0=0; i0<chunkSize; i0++) {
628              if (totalEle >= numEle) break; /* End inner loop */
629              fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
630              for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
631              fscanf(fileHandle_p, "\n");
632              totalEle++;
633              chunkEle++;
634            }
635    #ifdef PASO_MPI
636            /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
637            if (nextCPU < mpi_info->size) {
638                  int mpi_error;
639              tempInts[chunkSize*(2+numNodes)] = chunkEle;
640              mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81724, mpi_info->comm);
641              if ( mpi_error != MPI_SUCCESS ) {
642                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for ContactElements failed");
643                    return NULL;
644              }
645            }
646    #endif
647            nextCPU++;
648            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
649            if (nextCPU > mpi_info->size) break; /* End infinite loop */
650          } /* Infinite loop */
651        }   /* End master */
652        else {  /* Worker */
653    #ifdef PASO_MPI
654          /* Each worker receives one message */
655          MPI_Status status;
656          int mpi_error;
657          mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81724, mpi_info->comm, &status);
658          if ( mpi_error != MPI_SUCCESS ) {
659                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for ContactElements failed");
660                return NULL;
661          }
662          chunkEle = tempInts[chunkSize*(2+numNodes)];
663  #endif  #endif
664        }   /* Worker */
665    
666  #if 0 /* comment out the rest of the un-implemented crap for now */      /* Copy Element data from tempInts to mesh_p */
667          /* get the face elements */      Finley_ElementFile_allocTable(mesh_p->ContactElements, chunkEle);
668          if (Finley_noError()) {          mesh_p->ContactElements->minColor=0;
669               fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);          mesh_p->ContactElements->maxColor=chunkEle-1;
              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);  
              } else {  
                 mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);  
                 if (Finley_noError()) {  
                    Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle);  
                    if (Finley_noError()) {  
                       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 */  
670          if (Finley_noError()) {          if (Finley_noError()) {
671               fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);            #pragma omp parallel for private (i0, i1)
672               contactTypeID=Finley_RefElement_getTypeId(element_type);        for (i0=0; i0<chunkEle; i0++) {
673               if (contactTypeID==NoType) {          mesh_p->ContactElements->Id[i0] = tempInts[i0*(2+numNodes)+0];
674                 sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type);          mesh_p->ContactElements->Tag[i0]    = tempInts[i0*(2+numNodes)+1];
675                 Finley_setError(VALUE_ERROR,error_msg);              mesh_p->ContactElements->Color[i0] = i0;
676               } else {          for (i1 = 0; i1 < numNodes; i1++) {
677                 mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);            mesh_p->ContactElements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
678                 if (Finley_noError()) {          }
679                     Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);            }
                    if (Finley_noError()) {  
                       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 */  
                   }  
                }  
              }  
680          }          }
681          /* get the nodal element */  
682        TMPMEMFREE(tempInts);
683          } /* end of Read the contact element data */
684    
685            /* read nodal elements */
686    
687        /* Read the element typeID */
688          if (Finley_noError()) {          if (Finley_noError()) {
689               fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);        if (mpi_info->rank == 0) {
690               pointTypeID=Finley_RefElement_getTypeId(element_type);              fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
691               if (pointTypeID==NoType) {              typeID=Finley_RefElement_getTypeId(element_type);
692                 sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for points",element_type);        }
693                 Finley_setError(VALUE_ERROR,error_msg);  #ifdef PASO_MPI
694               }        if (mpi_info->size > 1) {
695               mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);          int temp1[2], mpi_error;
696               if (Finley_noError()) {          temp1[0] = (int) typeID;
697                  Finley_ElementFile_allocTable(mesh_p->Points, numEle);          temp1[1] = numEle;
698                  if (Finley_noError()) {          mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
699                     mesh_p->Points->minColor=0;          if (mpi_error != MPI_SUCCESS) {
700                     mesh_p->Points->maxColor=numEle-1;            Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
701                     for (i0 = 0; i0 < numEle; i0++) {            return NULL;
702                       fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);          }
703                       mesh_p->Points->Color[i0]=i0;          typeID = (ElementTypeId) temp1[0];
704                       for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {          numEle = temp1[1];
705                           fscanf(fileHandle_p, " %d",        }
706                              &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);  #endif
707                       }  /* for i1 */            if (typeID==NoType) {
708                       fscanf(fileHandle_p, "\n");              sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
709                     } /* for i0 */              Finley_setError(VALUE_ERROR, error_msg);
710                  }            }
711               }      }
712          }  
713          /* get the name tags */        /* Allocate the ElementFile */
714          mesh_p->Points=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
715          numNodes = mesh_p->Points->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
716    
717          /* Read the nodal element data */
718          if (Finley_noError()) {
719        int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
720        int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
721        /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
722        if (mpi_info->rank == 0) {  /* Master */
723          for (;;) {            /* Infinite loop */
724            for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
725            chunkEle = 0;
726            for (i0=0; i0<chunkSize; i0++) {
727              if (totalEle >= numEle) break; /* End inner loop */
728              fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
729              for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
730              fscanf(fileHandle_p, "\n");
731              totalEle++;
732              chunkEle++;
733            }
734    #ifdef PASO_MPI
735            /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
736            if (nextCPU < mpi_info->size) {
737                  int mpi_error;
738              tempInts[chunkSize*(2+numNodes)] = chunkEle;
739              mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81725, mpi_info->comm);
740              if ( mpi_error != MPI_SUCCESS ) {
741                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Points failed");
742                    return NULL;
743              }
744            }
745    #endif
746            nextCPU++;
747            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
748            if (nextCPU > mpi_info->size) break; /* End infinite loop */
749          } /* Infinite loop */
750        }   /* End master */
751        else {  /* Worker */
752    #ifdef PASO_MPI
753          /* Each worker receives one message */
754          MPI_Status status;
755          int mpi_error;
756          mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81725, mpi_info->comm, &status);
757          if ( mpi_error != MPI_SUCCESS ) {
758                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for Points failed");
759                return NULL;
760          }
761          chunkEle = tempInts[chunkSize*(2+numNodes)];
762    #endif
763        }   /* Worker */
764    
765        /* Copy Element data from tempInts to mesh_p */
766        Finley_ElementFile_allocTable(mesh_p->Points, chunkEle);
767            mesh_p->Points->minColor=0;
768            mesh_p->Points->maxColor=chunkEle-1;
769          if (Finley_noError()) {          if (Finley_noError()) {
770             if (feof(fileHandle_p) == 0) {            #pragma omp parallel for private (i0, i1)
771                fscanf(fileHandle_p, "%s\n", name);        for (i0=0; i0<chunkEle; i0++) {
772                while (feof(fileHandle_p) == 0) {          mesh_p->Points->Id[i0]  = tempInts[i0*(2+numNodes)+0];
773                     fscanf(fileHandle_p, "%s %d\n", name, &tag_key);          mesh_p->Points->Tag[i0] = tempInts[i0*(2+numNodes)+1];
774                     Finley_Mesh_addTagMap(mesh_p,name,tag_key);              mesh_p->Points->Color[i0] = i0;
775                }          for (i1 = 0; i1 < numNodes; i1++) {
776             }            mesh_p->Points->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
777            }
778              }
779          }          }
780  #endif /* comment out the rest of the un-implemented crap for now */  
781        TMPMEMFREE(tempInts);
782          } /* end of Read the nodal element data */
783    
784          /* ksteube TODO: read tags */
785    
786       }       }
787    
788       /* close file */       /* close file */
789       if (mpi_info->rank == 0) fclose(fileHandle_p);       if (mpi_info->rank == 0) fclose(fileHandle_p);
790    
791       /*   resolve id's : */       /*   resolve id's : */
792       /* rearrange elements: */       /* rearrange elements: */
793    
794         /* return mesh_p; */ /* ksteube temp return for debugging */
795    
796       if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);       if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
797       if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);       if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
      return mesh_p; /* ksteube temp return for debugging */  
798    
799       /* that's it */       /* that's it */
800       #ifdef Finley_TRACE       #ifdef Finley_TRACE

Legend:
Removed from v.1628  
changed lines
  Added in v.1732

  ViewVC Help
Powered by ViewVC 1.1.26