/[escript]/branches/arrayview_from_1695_trunk/finley/src/Mesh_read.c
ViewVC logotype

Diff of /branches/arrayview_from_1695_trunk/finley/src/Mesh_read.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1780 by jfenwick, Mon Aug 11 03:33:40 2008 UTC revision 1781 by jfenwick, Thu Sep 11 05:03:14 2008 UTC
# Line 19  Line 19 
19    
20  /**************************************************************/  /**************************************************************/
21    
22    #include <ctype.h>
23  #include "Mesh.h"  #include "Mesh.h"
24    
25  /**************************************************************/  /**************************************************************/
# Line 106  Finley_Mesh* Finley_Mesh_read(char* fnam Line 107  Finley_Mesh* Finley_Mesh_read(char* fnam
107                      for (i0 = 0; i0 < numEle; i0++) {                      for (i0 = 0; i0 < numEle; i0++) {
108                        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]);
109                        mesh_p->Elements->Color[i0]=i0;                        mesh_p->Elements->Color[i0]=i0;
110                          mesh_p->Elements->Owner[i0]=0;
111                        for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {                        for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {
112                             fscanf(fileHandle_p, " %d",                             fscanf(fileHandle_p, " %d",
113                                &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 135  Finley_Mesh* Finley_Mesh_read(char* fnam
135                        for (i0 = 0; i0 < numEle; i0++) {                        for (i0 = 0; i0 < numEle; i0++) {
136                          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]);
137                          mesh_p->FaceElements->Color[i0]=i0;                          mesh_p->FaceElements->Color[i0]=i0;
138                            mesh_p->FaceElements->Owner[i0]=0;
139                          for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {                          for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {
140                               fscanf(fileHandle_p, " %d",                               fscanf(fileHandle_p, " %d",
141                                  &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 163  Finley_Mesh* Finley_Mesh_read(char* fnam
163                        for (i0 = 0; i0 < numEle; i0++) {                        for (i0 = 0; i0 < numEle; i0++) {
164                          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]);
165                          mesh_p->ContactElements->Color[i0]=i0;                          mesh_p->ContactElements->Color[i0]=i0;
166                            mesh_p->ContactElements->Owner[i0]=0;
167                          for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {                          for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {
168                              fscanf(fileHandle_p, " %d",                              fscanf(fileHandle_p, " %d",
169                                 &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 191  Finley_Mesh* Finley_Mesh_read(char* fnam
191                     for (i0 = 0; i0 < numEle; i0++) {                     for (i0 = 0; i0 < numEle; i0++) {
192                       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]);
193                       mesh_p->Points->Color[i0]=i0;                       mesh_p->Points->Color[i0]=i0;
194                         mesh_p->Points->Owner[i0]=0;
195                       for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {                       for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {
196                           fscanf(fileHandle_p, " %d",                           fscanf(fileHandle_p, " %d",
197                              &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 241  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 246  Finley_Mesh* Finley_Mesh_read_MPI(char*
246    char error_msg[LenErrorMsg_MAX];    char error_msg[LenErrorMsg_MAX];
247    double time0=Finley_timer();    double time0=Finley_timer();
248    FILE *fileHandle_p = NULL;    FILE *fileHandle_p = NULL;
249    ElementTypeId typeID;    ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;
250      Finley_TagMap* tag_map;
 #if 0 /* comment out the rest of the un-implemented crap for now */  
       /* See below */  
   ElementTypeId faceTypeID, contactTypeID, pointTypeID;  
251    index_t tag_key;    index_t tag_key;
 #endif  
252    
253    Finley_resetError();    Finley_resetError();
254    
# Line 270  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 271  Finley_Mesh* Finley_Mesh_read_MPI(char*
271    }    }
272    
273  #ifdef PASO_MPI  #ifdef PASO_MPI
274    /* MPI Broadcast numDim, numNodes, name */    /* MPI Broadcast numDim, numNodes, name if there are multiple MPI procs*/
275    if (mpi_info->size > 0) {    if (mpi_info->size > 1) {
276      int temp1[3], error_code;      int temp1[3], error_code;
277      temp1[0] = numDim;      temp1[0] = numDim;
278      temp1[1] = numNodes;      temp1[1] = numNodes;
# Line 294  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 295  Finley_Mesh* Finley_Mesh_read_MPI(char*
295       /* allocate mesh */       /* allocate mesh */
296       mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);       mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
297       if (Finley_noError()) {       if (Finley_noError()) {
298        /* Each CPU will get at most chunkSize nodes so the message has to be sufficiently large */
299      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;
300      int *tempInts = TMPMEMALLOC(numNodes*3+1, index_t);      int *tempInts = TMPMEMALLOC(chunkSize*3+1, index_t);        /* Stores the integer message data */
301      double *tempCoords = TMPMEMALLOC(numNodes*numDim, double);      double *tempCoords = TMPMEMALLOC(chunkSize*numDim, double); /* Stores the double message data */
302    
303      /*      /*
304        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
305        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
306        First chunk sent to CPU 1, second to CPU 2, ...        First chunk sent to CPU 1, second to CPU 2, ...
307        Last chunk stays on CPU 0 (the master)        Last chunk stays on CPU 0 (the master)
# Line 308  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 310  Finley_Mesh* Finley_Mesh_read_MPI(char*
310    
311      if (mpi_info->rank == 0) {  /* Master */      if (mpi_info->rank == 0) {  /* Master */
312        for (;;) {            /* Infinite loop */        for (;;) {            /* Infinite loop */
313            for (i0=0; i0<chunkSize*3+1; i0++) tempInts[i0] = -1;
314            for (i0=0; i0<chunkSize*numDim; i0++) tempCoords[i0] = -1.0;
315          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;  
316          for (i1=0; i1<chunkSize; i1++) {          for (i1=0; i1<chunkSize; i1++) {
317            if (totalNodes >= numNodes) break;            if (totalNodes >= numNodes) break;    /* End of inner loop */
318                if (1 == numDim)                if (1 == numDim)
319          fscanf(fileHandle_p, "%d %d %d %le\n",          fscanf(fileHandle_p, "%d %d %d %le\n",
320            &tempInts[0+i1], &tempInts[numNodes+i1], &tempInts[numNodes*2+i1],            &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
321            &tempCoords[i1*numDim+0]);            &tempCoords[i1*numDim+0]);
322                if (2 == numDim)                if (2 == numDim)
323          fscanf(fileHandle_p, "%d %d %d %le %le\n",          fscanf(fileHandle_p, "%d %d %d %le %le\n",
324            &tempInts[0+i1], &tempInts[numNodes+i1], &tempInts[numNodes*2+i1],            &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
325            &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1]);            &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1]);
326                if (3 == numDim)                if (3 == numDim)
327          fscanf(fileHandle_p, "%d %d %d %le %le %le\n",          fscanf(fileHandle_p, "%d %d %d %le %le %le\n",
328            &tempInts[0+i1], &tempInts[numNodes+i1], &tempInts[numNodes*2+i1],            &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
329            &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1], &tempCoords[i1*numDim+2]);            &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1], &tempCoords[i1*numDim+2]);
330            totalNodes++;            totalNodes++; /* When do we quit the infinite loop? */
331            chunkNodes++;            chunkNodes++; /* How many nodes do we actually have in this chunk? It may be smaller than chunkSize. */
332            }
333            if (chunkNodes > chunkSize) {
334                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: error reading chunks of mesh, data too large for message size");
335                  return NULL;
336          }          }
         /* 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) {  
337  #ifdef PASO_MPI  #ifdef PASO_MPI
338            /* Eventually we'll send chunkSize nodes to each CPU numbered 1 ... mpi_info->size-1, here goes one of them */
339            if (nextCPU < mpi_info->size) {
340                int mpi_error;                int mpi_error;
341              tempInts[chunkSize*3] = chunkNodes;   /* The message has one more int to send chunkNodes */
342            tempInts[numNodes*3] = chunkNodes;            mpi_error = MPI_Send(tempInts, chunkSize*3+1, MPI_INT, nextCPU, 81720, mpi_info->comm);
           /* ksteube The size of this message can and should be brought down to chunkNodes*3+1, must re-org tempInts */  
           mpi_error = MPI_Send(tempInts, numNodes*3+1, MPI_INT, nextCPU, 81720, mpi_info->comm);  
343            if ( mpi_error != MPI_SUCCESS ) {            if ( mpi_error != MPI_SUCCESS ) {
344                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts failed");                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts failed");
345                  return NULL;                  return NULL;
346            }            }
347            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);
348            if ( mpi_error != MPI_SUCCESS ) {            if ( mpi_error != MPI_SUCCESS ) {
349                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempCoords failed");                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempCoords failed");
350                  return NULL;                  return NULL;
351            }            }
 #endif  
           nextCPU++;  
352          }          }
353          if (totalNodes >= numNodes) break;  #endif
354            nextCPU++;
355            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
356            if (nextCPU > mpi_info->size) break; /* End infinite loop */
357        } /* Infinite loop */        } /* Infinite loop */
358      }   /* End master */      }   /* End master */
359      else {  /* Worker */      else {  /* Worker */
# Line 356  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 361  Finley_Mesh* Finley_Mesh_read_MPI(char*
361        /* Each worker receives two messages */        /* Each worker receives two messages */
362        MPI_Status status;        MPI_Status status;
363            int mpi_error;            int mpi_error;
364        mpi_error = MPI_Recv(tempInts, numNodes*3+1, MPI_INT, 0, 81720, mpi_info->comm, &status);        mpi_error = MPI_Recv(tempInts, chunkSize*3+1, MPI_INT, 0, 81720, mpi_info->comm, &status);
365        if ( mpi_error != MPI_SUCCESS ) {        if ( mpi_error != MPI_SUCCESS ) {
366              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts failed");              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts failed");
367              return NULL;              return NULL;
368        }        }
369        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);
370        if ( mpi_error != MPI_SUCCESS ) {        if ( mpi_error != MPI_SUCCESS ) {
371              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempCoords failed");              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempCoords failed");
372              return NULL;              return NULL;
373        }        }
374        chunkNodes = tempInts[numNodes*3];        chunkNodes = tempInts[chunkSize*3];   /* How many nodes are in this workers chunk? */
375  #endif  #endif
376      }   /* Worker */      }   /* Worker */
377    
 #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);  
378      /* Copy node data from tempMem to mesh_p */      /* Copy node data from tempMem to mesh_p */
379          Finley_NodeFile_allocTable(mesh_p->Nodes, chunkNodes);          Finley_NodeFile_allocTable(mesh_p->Nodes, chunkNodes);
380          if (Finley_noError()) {          if (Finley_noError()) {
381        for (i0=0; i0<chunkNodes; i0++) {        for (i0=0; i0<chunkNodes; i0++) {
382          mesh_p->Nodes->Id[i0]               = tempInts[0+i0];          mesh_p->Nodes->Id[i0]               = tempInts[0+i0];
383          mesh_p->Nodes->globalDegreesOfFreedom[i0]       = tempInts[numNodes+i0];          mesh_p->Nodes->globalDegreesOfFreedom[i0]       = tempInts[chunkSize+i0];
384          mesh_p->Nodes->Tag[i0]              = tempInts[numNodes*2+i0];          mesh_p->Nodes->Tag[i0]              = tempInts[chunkSize*2+i0];
385          for (i1=0; i1<numDim; i1++) {          for (i1=0; i1<numDim; i1++) {
386            mesh_p->Nodes->Coordinates[INDEX2(i1,i0,numDim)]  = tempCoords[i0*numDim+i1];            mesh_p->Nodes->Coordinates[INDEX2(i1,i0,numDim)]  = tempCoords[i0*numDim+i1];
387          }          }
# Line 410  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 400  Finley_Mesh* Finley_Mesh_read_MPI(char*
400              typeID=Finley_RefElement_getTypeId(element_type);              typeID=Finley_RefElement_getTypeId(element_type);
401        }        }
402  #ifdef PASO_MPI  #ifdef PASO_MPI
403        if (mpi_info->size > 0) {        if (mpi_info->size > 1) {
404          int temp1[3], mpi_error;          int temp1[2], mpi_error;
405          temp1[0] = (int) typeID;          temp1[0] = (int) typeID;
406          temp1[1] = numEle;          temp1[1] = numEle;
407          mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);          mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
# Line 429  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 419  Finley_Mesh* Finley_Mesh_read_MPI(char*
419            }            }
420      }      }
421    
422        /* Read the element data */        /* Allocate the ElementFile */
423        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);
424        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 */
425    
426          /* Read the element data */
427        if (Finley_noError()) {        if (Finley_noError()) {
428      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;
429      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) */
430      if (numEle % mpi_info->size != 0) chunkSize++; /* Remainder from numEle / mpi_info->size will be spread out one-per-CPU */      /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
431      if (mpi_info->rank == 0) {  /* Master */      if (mpi_info->rank == 0) {  /* Master */
432        for (;;) {            /* Infinite loop */        for (;;) {            /* Infinite loop */
433            for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
434          chunkEle = 0;          chunkEle = 0;
         for (i0=0; i0<numEle*(2+numNodes)+1; i0++) tempInts[i0] = -1;  
435          for (i0=0; i0<chunkSize; i0++) {          for (i0=0; i0<chunkSize; i0++) {
436            if (totalEle >= numEle) break; /* End infinite loop */            if (totalEle >= numEle) break; /* End inner loop */
437            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]);
438            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]);
439            fscanf(fileHandle_p, "\n");            fscanf(fileHandle_p, "\n");
440            totalEle++;            totalEle++;
441            chunkEle++;            chunkEle++;
442          }          }
         /* Eventually we'll send chunk of nodes to each CPU except 0 itself, here goes one of them */  
         if (nextCPU < mpi_info->size) {  
443  #ifdef PASO_MPI  #ifdef PASO_MPI
444            /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
445            if (nextCPU < mpi_info->size) {
446                int mpi_error;                int mpi_error;
447              tempInts[chunkSize*(2+numNodes)] = chunkEle;
448            tempInts[numEle*(2+numNodes)] = chunkEle;            mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81722, mpi_info->comm);
 printf("ksteube CPU=%d/%d send to %d\n", mpi_info->rank, mpi_info->size, nextCPU);  
           mpi_error = MPI_Send(tempInts, numEle*(2+numNodes)+1, MPI_INT, nextCPU, 81722, mpi_info->comm);  
449            if ( mpi_error != MPI_SUCCESS ) {            if ( mpi_error != MPI_SUCCESS ) {
450                  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");
451                  return NULL;                  return NULL;
452            }            }
 #endif  
           nextCPU++;  
453          }          }
454          if (totalEle >= numEle) break; /* End infinite loop */  #endif
455            nextCPU++;
456            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
457            if (nextCPU > mpi_info->size) break; /* End infinite loop */
458        } /* Infinite loop */        } /* Infinite loop */
459      }   /* End master */      }   /* End master */
460      else {  /* Worker */      else {  /* Worker */
461  #ifdef PASO_MPI  #ifdef PASO_MPI
462        /* Each worker receives two messages */        /* Each worker receives one message */
463        MPI_Status status;        MPI_Status status;
464        int mpi_error;        int mpi_error;
465  printf("ksteube CPU=%d/%d recv on %d\n", mpi_info->rank, mpi_info->size, mpi_info->rank);        mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81722, mpi_info->comm, &status);
       mpi_error = MPI_Recv(tempInts, numEle*(2+numNodes)+1, MPI_INT, 0, 81722, mpi_info->comm, &status);  
466        if ( mpi_error != MPI_SUCCESS ) {        if ( mpi_error != MPI_SUCCESS ) {
467              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");
468              return NULL;              return NULL;
469        }        }
470        chunkEle = tempInts[numEle*(2+numNodes)];        chunkEle = tempInts[chunkSize*(2+numNodes)];
471  #endif  #endif
472      }   /* 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  
473    
474      /* Copy Element data from tempInts to mesh_p */      /* Copy Element data from tempInts to mesh_p */
475      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 480  printf("ksteube CPU=%d/%d recv on %d\n",
480        for (i0=0; i0<chunkEle; i0++) {        for (i0=0; i0<chunkEle; i0++) {
481          mesh_p->Elements->Id[i0]    = tempInts[i0*(2+numNodes)+0];          mesh_p->Elements->Id[i0]    = tempInts[i0*(2+numNodes)+0];
482          mesh_p->Elements->Tag[i0]   = tempInts[i0*(2+numNodes)+1];          mesh_p->Elements->Tag[i0]   = tempInts[i0*(2+numNodes)+1];
483                mesh_p->Elements->Owner[i0]  =mpi_info->rank;
484                mesh_p->Elements->Color[i0] = i0;
485          for (i1 = 0; i1 < numNodes; i1++) {          for (i1 = 0; i1 < numNodes; i1++) {
486            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];
487          }          }
# Line 506  printf("ksteube CPU=%d/%d recv on %d\n", Line 489  printf("ksteube CPU=%d/%d recv on %d\n",
489          }          }
490    
491      TMPMEMFREE(tempInts);      TMPMEMFREE(tempInts);
492        }        } /* end of Read the element data */
493    
494            /* read face elements */
495    
496  #if 0 /* this is the original code for reading elements */      /* Read the element typeID */
         /* read elements */  
497          if (Finley_noError()) {          if (Finley_noError()) {
498          if (mpi_info->rank == 0) {
499                fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
500                typeID=Finley_RefElement_getTypeId(element_type);
501          }
502    #ifdef PASO_MPI
503          if (mpi_info->size > 1) {
504            int temp1[2], mpi_error;
505            temp1[0] = (int) typeID;
506            temp1[1] = numEle;
507            mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
508            if (mpi_error != MPI_SUCCESS) {
509              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
510              return NULL;
511            }
512            typeID = (ElementTypeId) temp1[0];
513            numEle = temp1[1];
514          }
515    #endif
516              if (typeID==NoType) {
517                sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
518                Finley_setError(VALUE_ERROR, error_msg);
519              }
520        }
521    
522             fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);        /* Allocate the ElementFile */
523             typeID=Finley_RefElement_getTypeId(element_type);        mesh_p->FaceElements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
524             if (typeID==NoType) {        numNodes = mesh_p->FaceElements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
525               sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s",element_type);  
526               Finley_setError(VALUE_ERROR,error_msg);        /* Read the face element data */
527             } else {        if (Finley_noError()) {
528               /* read the elements */      int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
529               mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);      int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
530               if (Finley_noError()) {      /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
531                   Finley_ElementFile_allocTable(mesh_p->Elements, numEle);      if (mpi_info->rank == 0) {  /* Master */
532                   mesh_p->Elements->minColor=0;        for (;;) {            /* Infinite loop */
533                   mesh_p->Elements->maxColor=numEle-1;          for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
534                   if (Finley_noError()) {          chunkEle = 0;
535                      for (i0 = 0; i0 < numEle; i0++) {          for (i0=0; i0<chunkSize; i0++) {
536                        fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]);            if (totalEle >= numEle) break; /* End inner loop */
537                        mesh_p->Elements->Color[i0]=i0;            fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
538                        for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {            for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
539                             fscanf(fileHandle_p, " %d",            fscanf(fileHandle_p, "\n");
540                                &mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);            totalEle++;
541                        } /* for i1 */            chunkEle++;
542                        fscanf(fileHandle_p, "\n");          }
543                      } /* for i0 */  #ifdef PASO_MPI
544                   }          /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
545               }          if (nextCPU < mpi_info->size) {
546                  int mpi_error;
547              tempInts[chunkSize*(2+numNodes)] = chunkEle;
548              mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81723, mpi_info->comm);
549              if ( mpi_error != MPI_SUCCESS ) {
550                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for FaceElements failed");
551                    return NULL;
552              }
553            }
554    #endif
555            nextCPU++;
556            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
557            if (nextCPU > mpi_info->size) break; /* End infinite loop */
558          } /* Infinite loop */
559        }   /* End master */
560        else {  /* Worker */
561    #ifdef PASO_MPI
562          /* Each worker receives one message */
563          MPI_Status status;
564          int mpi_error;
565          mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81723, mpi_info->comm, &status);
566          if ( mpi_error != MPI_SUCCESS ) {
567                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for FaceElements failed");
568                return NULL;
569          }
570          chunkEle = tempInts[chunkSize*(2+numNodes)];
571    #endif
572        }   /* Worker */
573    
574        /* Copy Element data from tempInts to mesh_p */
575        Finley_ElementFile_allocTable(mesh_p->FaceElements, chunkEle);
576            mesh_p->FaceElements->minColor=0;
577            mesh_p->FaceElements->maxColor=chunkEle-1;
578            if (Finley_noError()) {
579              #pragma omp parallel for private (i0, i1)
580          for (i0=0; i0<chunkEle; i0++) {
581            mesh_p->FaceElements->Id[i0]    = tempInts[i0*(2+numNodes)+0];
582            mesh_p->FaceElements->Tag[i0]   = tempInts[i0*(2+numNodes)+1];
583                mesh_p->FaceElements->Owner[i0]  =mpi_info->rank;
584                mesh_p->FaceElements->Color[i0] = i0;
585            for (i1 = 0; i1 < numNodes; i1++) {
586              mesh_p->FaceElements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
587            }
588            }            }
589          }          }
590    
591        TMPMEMFREE(tempInts);
592          } /* end of Read the face element data */
593    
594            /* read contact elements */
595    
596        /* Read the element typeID */
597            if (Finley_noError()) {
598          if (mpi_info->rank == 0) {
599                fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
600                typeID=Finley_RefElement_getTypeId(element_type);
601          }
602    #ifdef PASO_MPI
603          if (mpi_info->size > 1) {
604            int temp1[2], mpi_error;
605            temp1[0] = (int) typeID;
606            temp1[1] = numEle;
607            mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
608            if (mpi_error != MPI_SUCCESS) {
609              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
610              return NULL;
611            }
612            typeID = (ElementTypeId) temp1[0];
613            numEle = temp1[1];
614          }
615  #endif  #endif
616              if (typeID==NoType) {
617                sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
618                Finley_setError(VALUE_ERROR, error_msg);
619              }
620        }
621    
622  printf("ksteube CPU=%d/%d Element typeID=%d\n", mpi_info->rank, mpi_info->size, typeID);        /* Allocate the ElementFile */
623          mesh_p->ContactElements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
624          numNodes = mesh_p->ContactElements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
625    
626          /* Read the contact element data */
627          if (Finley_noError()) {
628        int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
629        int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
630        /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
631        if (mpi_info->rank == 0) {  /* Master */
632          for (;;) {            /* Infinite loop */
633            for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
634            chunkEle = 0;
635            for (i0=0; i0<chunkSize; i0++) {
636              if (totalEle >= numEle) break; /* End inner loop */
637              fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
638              for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
639              fscanf(fileHandle_p, "\n");
640              totalEle++;
641              chunkEle++;
642            }
643    #ifdef PASO_MPI
644            /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
645            if (nextCPU < mpi_info->size) {
646                  int mpi_error;
647              tempInts[chunkSize*(2+numNodes)] = chunkEle;
648              mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81724, mpi_info->comm);
649              if ( mpi_error != MPI_SUCCESS ) {
650                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for ContactElements failed");
651                    return NULL;
652              }
653            }
654    #endif
655            nextCPU++;
656            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
657            if (nextCPU > mpi_info->size) break; /* End infinite loop */
658          } /* Infinite loop */
659        }   /* End master */
660        else {  /* Worker */
661    #ifdef PASO_MPI
662          /* Each worker receives one message */
663          MPI_Status status;
664          int mpi_error;
665          mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81724, mpi_info->comm, &status);
666          if ( mpi_error != MPI_SUCCESS ) {
667                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for ContactElements failed");
668                return NULL;
669          }
670          chunkEle = tempInts[chunkSize*(2+numNodes)];
671    #endif
672        }   /* Worker */
673    
674  #if 1      /* Copy Element data from tempInts to mesh_p */
675        Finley_ElementFile_allocTable(mesh_p->ContactElements, chunkEle);
676            mesh_p->ContactElements->minColor=0;
677            mesh_p->ContactElements->maxColor=chunkEle-1;
678            if (Finley_noError()) {
679              #pragma omp parallel for private (i0, i1)
680          for (i0=0; i0<chunkEle; i0++) {
681            mesh_p->ContactElements->Id[i0] = tempInts[i0*(2+numNodes)+0];
682            mesh_p->ContactElements->Tag[i0]    = tempInts[i0*(2+numNodes)+1];
683                mesh_p->ContactElements->Owner[i0]  =mpi_info->rank;
684                mesh_p->ContactElements->Color[i0] = i0;
685            for (i1 = 0; i1 < numNodes; i1++) {
686              mesh_p->ContactElements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
687            }
688              }
689            }
690    
691    /* Define other structures to keep mesh_write from crashing */      TMPMEMFREE(tempInts);
692    /* Change the typeid from NoType later */        } /* end of Read the contact element data */
693    
694    mesh_p->FaceElements=Finley_ElementFile_alloc(NoType, mesh_p->order, mesh_p->reduced_order, mpi_info);          /* read nodal elements */
   Finley_ElementFile_allocTable(mesh_p->FaceElements, 0);  
695    
696    mesh_p->ContactElements=Finley_ElementFile_alloc(NoType, mesh_p->order, mesh_p->reduced_order, mpi_info);      /* Read the element typeID */
697    Finley_ElementFile_allocTable(mesh_p->ContactElements, 0);          if (Finley_noError()) {
698          if (mpi_info->rank == 0) {
699                fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
700                typeID=Finley_RefElement_getTypeId(element_type);
701          }
702    #ifdef PASO_MPI
703          if (mpi_info->size > 1) {
704            int temp1[2], mpi_error;
705            temp1[0] = (int) typeID;
706            temp1[1] = numEle;
707            mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
708            if (mpi_error != MPI_SUCCESS) {
709              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
710              return NULL;
711            }
712            typeID = (ElementTypeId) temp1[0];
713            numEle = temp1[1];
714          }
715    #endif
716              if (typeID==NoType) {
717                sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
718                Finley_setError(VALUE_ERROR, error_msg);
719              }
720        }
721    
722    mesh_p->Points=Finley_ElementFile_alloc(NoType, mesh_p->order, mesh_p->reduced_order, mpi_info);        /* Allocate the ElementFile */
723    Finley_ElementFile_allocTable(mesh_p->Points, 0);        mesh_p->Points=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
724          numNodes = mesh_p->Points->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
725    
726          /* Read the nodal element data */
727          if (Finley_noError()) {
728        int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
729        int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
730        /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
731        if (mpi_info->rank == 0) {  /* Master */
732          for (;;) {            /* Infinite loop */
733            for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
734            chunkEle = 0;
735            for (i0=0; i0<chunkSize; i0++) {
736              if (totalEle >= numEle) break; /* End inner loop */
737              fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
738              for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
739              fscanf(fileHandle_p, "\n");
740              totalEle++;
741              chunkEle++;
742            }
743    #ifdef PASO_MPI
744            /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
745            if (nextCPU < mpi_info->size) {
746                  int mpi_error;
747              tempInts[chunkSize*(2+numNodes)] = chunkEle;
748              mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81725, mpi_info->comm);
749              if ( mpi_error != MPI_SUCCESS ) {
750                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Points failed");
751                    return NULL;
752              }
753            }
754  #endif  #endif
755            nextCPU++;
756            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
757            if (nextCPU > mpi_info->size) break; /* End infinite loop */
758          } /* Infinite loop */
759        }   /* End master */
760        else {  /* Worker */
761    #ifdef PASO_MPI
762          /* Each worker receives one message */
763          MPI_Status status;
764          int mpi_error;
765          mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81725, mpi_info->comm, &status);
766          if ( mpi_error != MPI_SUCCESS ) {
767                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for Points failed");
768                return NULL;
769          }
770          chunkEle = tempInts[chunkSize*(2+numNodes)];
771    #endif
772        }   /* Worker */
773    
774  #if 0 /* comment out the rest of the un-implemented crap for now */      /* Copy Element data from tempInts to mesh_p */
775          /* get the face elements */      Finley_ElementFile_allocTable(mesh_p->Points, chunkEle);
776          if (Finley_noError()) {          mesh_p->Points->minColor=0;
777               fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);          mesh_p->Points->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 */  
778          if (Finley_noError()) {          if (Finley_noError()) {
779               fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);            #pragma omp parallel for private (i0, i1)
780               contactTypeID=Finley_RefElement_getTypeId(element_type);        for (i0=0; i0<chunkEle; i0++) {
781               if (contactTypeID==NoType) {          mesh_p->Points->Id[i0]  = tempInts[i0*(2+numNodes)+0];
782                 sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type);          mesh_p->Points->Tag[i0] = tempInts[i0*(2+numNodes)+1];
783                 Finley_setError(VALUE_ERROR,error_msg);              mesh_p->Points->Owner[i0]  =mpi_info->rank;
784               } else {              mesh_p->Points->Color[i0] = i0;
785                 mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);          for (i1 = 0; i1 < numNodes; i1++) {
786                 if (Finley_noError()) {            mesh_p->Points->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
787                     Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);          }
788                     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 */  
                   }  
                }  
              }  
789          }          }
790          /* get the nodal element */  
791          if (Finley_noError()) {      TMPMEMFREE(tempInts);
792               fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);        } /* end of Read the nodal element data */
793               pointTypeID=Finley_RefElement_getTypeId(element_type);  
794               if (pointTypeID==NoType) {        /* get the name tags */
795                 sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for points",element_type);        if (Finley_noError()) {
796                 Finley_setError(VALUE_ERROR,error_msg);          char *remainder, *ptr;
797               }          int tag_key, len, error_code;
798               mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);          long cur_pos, end_pos;
799               if (Finley_noError()) {          if (mpi_info->rank == 0) {  /* Master */
800                  Finley_ElementFile_allocTable(mesh_p->Points, numEle);        /* Read the word 'Tag' */
801                  if (Finley_noError()) {        if (! feof(fileHandle_p)) fscanf(fileHandle_p, "%s\n", name);
802                     mesh_p->Points->minColor=0;        /* Read rest of file in one chunk, after using seek to find length */
803                     mesh_p->Points->maxColor=numEle-1;            cur_pos = ftell(fileHandle_p);
804                     for (i0 = 0; i0 < numEle; i0++) {            fseek(fileHandle_p, 0L, SEEK_END);
805                       fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);            end_pos = ftell(fileHandle_p);
806                       mesh_p->Points->Color[i0]=i0;            fseek(fileHandle_p, (long)cur_pos, SEEK_SET);
807                       for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {        remainder = TMPMEMALLOC(end_pos-cur_pos+1, char);
808                           fscanf(fileHandle_p, " %d",        if (! feof(fileHandle_p)) fread(remainder, (size_t) end_pos-cur_pos, sizeof(char), fileHandle_p);
809                              &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);        remainder[end_pos-cur_pos] = 0;
810                       }  /* for i1 */        len = strlen(remainder);    
811                       fscanf(fileHandle_p, "\n");        while ((--len)>0 && isspace(remainder[len])) remainder[len]=0;
812                     } /* for i0 */        len = strlen(remainder);
                 }  
              }  
813          }          }
814          /* get the name tags */  #ifdef PASO_MPI
815          if (Finley_noError()) {          error_code = MPI_Bcast (&len, 1, MPI_INT,  0, mpi_info->comm);
816             if (feof(fileHandle_p) == 0) {          if (error_code != MPI_SUCCESS) {
817                fscanf(fileHandle_p, "%s\n", name);            Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of tag len failed");
818                while (feof(fileHandle_p) == 0) {            return NULL;
819                     fscanf(fileHandle_p, "%s %d\n", name, &tag_key);          }
820                     Finley_Mesh_addTagMap(mesh_p,name,tag_key);      if (mpi_info->rank != 0) {
821                }        remainder = TMPMEMALLOC(len+1, char);
822             }        remainder[0] = 0;
823        }
824            error_code = MPI_Bcast (remainder, len+1, MPI_CHAR,  0, mpi_info->comm);
825            if (error_code != MPI_SUCCESS) {
826              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of tags failed");
827              return NULL;
828          }          }
829  #endif /* comment out the rest of the un-implemented crap for now */  #endif
830        if (remainder[0]) {
831              ptr = remainder;
832              do {
833                sscanf(ptr, "%s %d\n", name, &tag_key);
834                if (*name) Finley_Mesh_addTagMap(mesh_p,name,tag_key);
835                ptr++;
836              } while(NULL != (ptr = strchr(ptr, '\n')) && *ptr);
837              TMPMEMFREE(remainder);
838        }
839          }
840    
841       }       }
842    
843       /* close file */       /* close file */
844       if (mpi_info->rank == 0) fclose(fileHandle_p);       if (mpi_info->rank == 0) fclose(fileHandle_p);
845    
# Line 660  printf("ksteube CPU=%d/%d Element typeID Line 848  printf("ksteube CPU=%d/%d Element typeID
848    
849       if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);       if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
850       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 */  
851    
852       /* that's it */       /* that's it */
853       #ifdef Finley_TRACE       #ifdef Finley_TRACE

Legend:
Removed from v.1780  
changed lines
  Added in v.1781

  ViewVC Help
Powered by ViewVC 1.1.26