/[escript]/trunk/finley/src/Mesh_saveVTK.c
ViewVC logotype

Diff of /trunk/finley/src/Mesh_saveVTK.c

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

revision 796 by dhawcroft, Mon Jul 31 06:16:17 2006 UTC revision 1028 by gross, Wed Mar 14 00:15:24 2007 UTC
# Line 43  Line 43 
43   out. The  routines are collectively called and will be called in the natural   out. The  routines are collectively called and will be called in the natural
44   ordering i.e 0 to maxProcs-1.   ordering i.e 0 to maxProcs-1.
45    
46     Notable Notables:
47     the struct localIndexCache: stores local domain indices for faster  reference
48  */  */
49    
50  #ifdef PASO_MPI  #ifdef PASO_MPI
# Line 60  Line 62 
62    
63  void Finley_Mesh_saveVTK_MPIO(const char * filename_p, Finley_Mesh *mesh_p, const dim_t num_data,char* *names_p, escriptDataC* *data_pp)  void Finley_Mesh_saveVTK_MPIO(const char * filename_p, Finley_Mesh *mesh_p, const dim_t num_data,char* *names_p, escriptDataC* *data_pp)
64  {  {
   double time0 = Paso_timer();  
65    int    numPoints,    int    numPoints,
66    numCells = -1,    numCells = -1,
67               myRank,comm,gsize,               myRank,comm,gsize,
# Line 68  void Finley_Mesh_saveVTK_MPIO(const char Line 69  void Finley_Mesh_saveVTK_MPIO(const char
69               nDim,               nDim,
70               shape;               shape;
71    size_t __n;    size_t __n;
   int* failSend;  
72    int i,j,k,m,n,count;    int i,j,k,m,n,count;
73    int numGlobalCells = 0;    int numGlobalCells = 0;
74    index_t  *nodesGlobal=NULL;   // used to get the connectivity  right for VTK    index_t  *nodesGlobal=NULL;   // used to get the connectivity  right for VTK
# Line 118  void Finley_Mesh_saveVTK_MPIO(const char Line 118  void Finley_Mesh_saveVTK_MPIO(const char
118    int numInternalNodes,    int numInternalNodes,
119    numLocalNodes,    numLocalNodes,
120    numBoundaryNodes,    numBoundaryNodes,
121    localDOF;    localDOF;  // equals to  (DOF of Internal Nodes) +  (DOF of Boundary Nodes) of local domain
122    
123    
124    nDim  = mesh_p->Nodes->numDim;    nDim  = mesh_p->Nodes->numDim;
# Line 139  void Finley_Mesh_saveVTK_MPIO(const char Line 139  void Finley_Mesh_saveVTK_MPIO(const char
139    infoHints = MPI_INFO_NULL;    infoHints = MPI_INFO_NULL;
140  #endif  #endif
141    
142    // Holds a local node/element index to their global arrays    // Holds a local node/element index into the global array
143    struct localIndexCache    struct localIndexCache
144    {    {
145      index_t *values;      index_t *values;
146      int size;      int size;
147    };    };
   typedef struct localIndexCache localIndexCache;  
148    
149    localIndexCache nodeCache,    struct localIndexCache nodeCache,
150    elementCache;            elementCache;
151    
152    // Collective Call    // Collective Call
153    MPI_File_open(mesh_p->Nodes->MPIInfo->comm, (char*)filename_p, amode,infoHints, &fh);    MPI_File_open(mesh_p->Nodes->MPIInfo->comm, (char*)filename_p, amode,infoHints, &fh);
# Line 453  void Finley_Mesh_saveVTK_MPIO(const char Line 452  void Finley_Mesh_saveVTK_MPIO(const char
452              "<UnstructuredGrid>\n" \              "<UnstructuredGrid>\n" \
453              "<Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n" \              "<Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n" \
454              "<Points>\n" \              "<Points>\n" \
455              "<DataArray NumberOfComponents=\"%d\" type=\"Float32\" format=\"ascii\">\n"              "<DataArray NumberOfComponents=\"%d\" type=\"Float64\" format=\"ascii\">\n"
456              ,numPoints,numGlobalCells,MAX(3,nDim));              ,numPoints,numGlobalCells,MAX(3,nDim));
457    
458    
# Line 465  void Finley_Mesh_saveVTK_MPIO(const char Line 464  void Finley_Mesh_saveVTK_MPIO(const char
464    
465    numLocalNodes=localDOF;    numLocalNodes=localDOF;
466    
467    //  values vary from 13-14 chars hence the strlen()    //  we will be writing values which vary from 13-15 chars hence the strlen()
468    char* largebuf = MEMALLOC( numLocalNodes*14*nDim + numLocalNodes*2 + 1 ,char);    char* largebuf = MEMALLOC( numLocalNodes*15*nDim + numLocalNodes*2 + 1 ,char);
469    largebuf[0] = '\0';    largebuf[0] = '\0';
470    char tmpbuf[14];    char tmpbuf[15];
471    int tsz=0;    int tsz=0;
472    int numNodesOutput=0;    int numNodesOutput=0;
473    index_t pos=0;    index_t pos=0;
# Line 479  void Finley_Mesh_saveVTK_MPIO(const char Line 478  void Finley_Mesh_saveVTK_MPIO(const char
478    nodeCache.values = MEMALLOC( numLocalNodes, index_t);    nodeCache.values = MEMALLOC( numLocalNodes, index_t);
479    index_t bc_pos = 0;    index_t bc_pos = 0;
480    
481    // Custom strcat:  avoids expensive strlen(3) call by  strcat(3)    // Custom string concat:  avoids expensive strlen(3) call by strcat(3)
482      // Note the implicit assumption on the variable "tsz"
483    int __len,__j;    int __len,__j;
484    char  *zero = "0.000000e+00 ";    char  *zero = "0.000000e+00 ";
485    char  *newline = "\n";    char  *newline = "\n";
486  #define __STRCAT(__buf,x)  \    
487    #define __STRCAT(dest,chunk,tsz)  \
488  {                  \  {                  \
489       __len = strlen(chunk); \
490     __j = -1;      \     __j = -1;      \
491     while(__j++ < __len)  \     while(__j++ < __len)  \
492      *(__buf+tsz+__j)=*(x+__j); \      *(dest+tsz+__j)=*(chunk+__j); \
493       tsz+=__len;              \
494  }  }
495      
496      // Loop over all nodes    
497    for (i = 0; i < mesh_p->Nodes->numNodes; i++)    for (i = 0; i < mesh_p->Nodes->numNodes; i++)
498    {    {
499      // This is the bit that will break for periodic BCs because it assumes that there is a one to one      // This is the bit that will break for periodic BCs because it assumes that there is a one to one
500      // correspondance between nodes and Degrees of freedom      // correspondance between nodes and Degrees of freedom
501        //TODO: handle periodic BC's
502      DOFNodes[mesh_p->Nodes->degreeOfFreedom[i]] = i;      DOFNodes[mesh_p->Nodes->degreeOfFreedom[i]] = i;
503    
504      /* local node ?*/      // Is this node local to the domain ?
505      if( mesh_p->Nodes->degreeOfFreedom[i] < localDOF )      if( mesh_p->Nodes->degreeOfFreedom[i] < localDOF )
506      {      {
507        for (j = 0; j < nDim; j++)        for (j = 0; j < nDim; j++)
508        {        {
509          sprintf(tmpbuf,"%e ", mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)] );          sprintf(tmpbuf,"%e ", mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)] );
510          __len = strlen(tmpbuf);          __STRCAT(largebuf,tmpbuf,tsz)
         __STRCAT(largebuf,tmpbuf)  
         tsz+=__len;  
511        }        }
512        for (k=0; k<3-nDim; k++)        for (k=0; k<3-nDim; k++)
513        {        {
514          __len = 13;          __STRCAT(largebuf,zero,tsz)
515          __STRCAT(largebuf,zero)        }
516          tsz+=__len;        __STRCAT(largebuf,newline,tsz)
       }  
       __len=1;  
       __STRCAT(largebuf,newline)  
       tsz += 1;  
517        nodeCache.values[numNodesOutput++]=i;        nodeCache.values[numNodesOutput++]=i;
518      }      }
519    }    }
# Line 568  void Finley_Mesh_saveVTK_MPIO(const char Line 567  void Finley_Mesh_saveVTK_MPIO(const char
567      {      {
568        Paso_CommBuffer_recv(mesh_p->Nodes->CommBuffer, dist->neighbours[n], sizeof(index_t));        Paso_CommBuffer_recv(mesh_p->Nodes->CommBuffer, dist->neighbours[n], sizeof(index_t));
569        Paso_CommBuffer_unpack(mesh_p->Nodes->CommBuffer, dist->neighbours[n], NULL, backwardBuffer, sizeof(index_t), 0 );        Paso_CommBuffer_unpack(mesh_p->Nodes->CommBuffer, dist->neighbours[n], NULL, backwardBuffer, sizeof(index_t), 0 );
570          /* TODO: voodoo to handle perdiodic  BC's */
571        for( i=0; i<dist->edges[n]->numBackward; i++ )        for( i=0; i<dist->edges[n]->numBackward; i++ )
572          nodesGlobal[DOFNodes[dist->edges[n]->indexBackward[i] ]] = backwardBuffer[i];          nodesGlobal[DOFNodes[dist->edges[n]->indexBackward[i] ]] = backwardBuffer[i];
573      }      }
574    }    }
575      
576    
577      
578    MEMFREE(vtxdist);    MEMFREE(vtxdist);
579    MEMFREE(DOFNodes);    MEMFREE(DOFNodes);
580    MEMFREE(backwardBuffer);    MEMFREE(backwardBuffer);
# Line 594  void Finley_Mesh_saveVTK_MPIO(const char Line 596  void Finley_Mesh_saveVTK_MPIO(const char
596    // Collective    // Collective
597    MPIO_DEBUG(" Writing Connectivity... ")    MPIO_DEBUG(" Writing Connectivity... ")
598    
   // TODO: Improve on upper bound , will fail for very large meshes!!  
599    size_t sz = numLocalCells*6*numVTKNodesPerElement + numLocalCells;    size_t sz = numLocalCells*6*numVTKNodesPerElement + numLocalCells;
600    largebuf = MEMALLOC(sz,char);    largebuf = MEMALLOC(sz,char);
601    largebuf[0] = '\0';    largebuf[0] = '\0';
# Line 612  void Finley_Mesh_saveVTK_MPIO(const char Line 613  void Finley_Mesh_saveVTK_MPIO(const char
613          for (j = 0; j < numVTKNodesPerElement; j++)          for (j = 0; j < numVTKNodesPerElement; j++)
614          {          {
615            sprintf(tmpbuf,"%d ",nodesGlobal[mesh_p->Nodes->toReduced[elements->Nodes[INDEX2(elements->ReferenceElement->Type->linearNodes[j], i, NN)]]]);            sprintf(tmpbuf,"%d ",nodesGlobal[mesh_p->Nodes->toReduced[elements->Nodes[INDEX2(elements->ReferenceElement->Type->linearNodes[j], i, NN)]]]);
616            __len=strlen(tmpbuf);            __STRCAT(largebuf,tmpbuf,tsz)
617            __STRCAT(largebuf,tmpbuf)          }
618            tsz+=__len;          __STRCAT(largebuf,newline,tsz)
         }  
         __len=1;  
         __STRCAT(largebuf,newline)  
         tsz+=1;  
   
619          elementCache.values[pos++]=i;          elementCache.values[pos++]=i;
620        }        }
621      }      }
622    }    }
623    else if (VTK_QUADRATIC_HEXAHEDRON==cellType)    else if (VTK_QUADRATIC_HEXAHEDRON==cellType)
624    {    {
625      char tmpbuf2[20*20];      char tmpbuf2[20*20*2];
626    
627      for (i = 0; i < numCells; i++)      for (i = 0; i < numCells; i++)
628      {      {
# Line 653  void Finley_Mesh_saveVTK_MPIO(const char Line 649  void Finley_Mesh_saveVTK_MPIO(const char
649                  nodesGlobal[elements->Nodes[INDEX2(13, i, NN)]],                  nodesGlobal[elements->Nodes[INDEX2(13, i, NN)]],
650                  nodesGlobal[elements->Nodes[INDEX2(14, i, NN)]],                  nodesGlobal[elements->Nodes[INDEX2(14, i, NN)]],
651                  nodesGlobal[elements->Nodes[INDEX2(15, i, NN)]]);                  nodesGlobal[elements->Nodes[INDEX2(15, i, NN)]]);
652          __len=strlen(tmpbuf2);          __STRCAT(largebuf,tmpbuf2,tsz)
         __STRCAT(largebuf,tmpbuf2)  
         tsz+=__len;  
653          elementCache.values[pos++]=i;          elementCache.values[pos++]=i;
654        }        }
655      }      }
# Line 665  void Finley_Mesh_saveVTK_MPIO(const char Line 659  void Finley_Mesh_saveVTK_MPIO(const char
659    
660      for (i = 0; i < numCells; i++)      for (i = 0; i < numCells; i++)
661      {      {
662        for (j = 0; j < numVTKNodesPerElement; j++)        if (elements->Id[i] >= elements->elementDistribution->vtxdist[i] &&  elements->Id[i] <= elements->elementDistribution->vtxdist[i+1] - 1 )
663        {        {
664          sprintf(tmpbuf,"%d ", nodesGlobal[elements->Nodes[INDEX2(elements->ReferenceElement->Type->geoNodes[j], i, NN)]]);          for (j = 0; j < numVTKNodesPerElement; j++)
665          __len=strlen(tmpbuf);          {
666          __STRCAT(largebuf,tmpbuf)            sprintf(tmpbuf,"%d ", nodesGlobal[elements->Nodes[INDEX2(elements->ReferenceElement->Type->geoNodes[j], i, NN)]]);
667          tsz+=__len;            __STRCAT(largebuf,tmpbuf,tsz)
668        }          }
669        __len=1;          __STRCAT(largebuf,newline,tsz)
670        __STRCAT(largebuf,newline)          elementCache.values[pos++]=i;
671        tsz+=1;        }
       elementCache.values[pos++]=i;  
672      }      }
673    }    }
674    else    else
675    {    {
   
676      for(i = 0;i  < numCells ; i++)      for(i = 0;i  < numCells ; i++)
677      {      {
       // is this element in domain of process with "myRank"  
   
678        if( elements->Id[i] >= elements->elementDistribution->vtxdist[myRank] && elements->Id[i] <= elements->elementDistribution->vtxdist[myRank+1]-1)        if( elements->Id[i] >= elements->elementDistribution->vtxdist[myRank] && elements->Id[i] <= elements->elementDistribution->vtxdist[myRank+1]-1)
679        {        {
680          for (j = 0; j < numVTKNodesPerElement; j++)          for (j = 0; j < numVTKNodesPerElement; j++)
681          {          {
682            sprintf(tmpbuf,"%d ", nodesGlobal[ elements->Nodes[INDEX2(j, i, NN) ] ] );            sprintf(tmpbuf,"%d ", nodesGlobal[ elements->Nodes[INDEX2(j, i, NN) ] ] );
683            __len=strlen(tmpbuf);            __STRCAT(largebuf,tmpbuf,tsz)
684            __STRCAT(largebuf,tmpbuf)          }
685            tsz+=__len;          __STRCAT(largebuf,newline,tsz)
         }  
         __len=1;  
         __STRCAT(largebuf,newline)  
         tsz+=1;  
686          elementCache.values[pos++]=i;          elementCache.values[pos++]=i;
687        }        }
688      }      }
689    }    }
690    
691    elementCache.size = pos;    elementCache.size = pos;
692    
693      largebuf[tsz] = '\0';
694    MPI_File_write_ordered(fh,largebuf,tsz, MPI_CHAR, &status);    MPI_File_write_ordered(fh,largebuf,tsz, MPI_CHAR, &status);
695    MEMFREE(largebuf);    MEMFREE(largebuf);
696    MPIO_DEBUG(" Done Writing Connectivity ")    MPIO_DEBUG(" Done Writing Connectivity ")
# Line 719  void Finley_Mesh_saveVTK_MPIO(const char Line 707  void Finley_Mesh_saveVTK_MPIO(const char
707    
708      int n = numVTKNodesPerElement;      int n = numVTKNodesPerElement;
709    
710        // allocate an upper bound on number of bytes needed  
711      int sz=0;      int sz=0;
712      int lg = log10(numGlobalCells * n) + 1;      int lg = log10(numGlobalCells * n) + 1;
713      sz += numGlobalCells*lg;      sz += numGlobalCells*lg;
714      sz += numGlobalCells;      sz += numGlobalCells;  
715      tsz = 0;      tsz = 0;
716    
717      char* largebuf = MEMALLOC(sz  + strlen(tag1) + strlen(tag2) + strlen(tag3) + strlen(tag4),char);      char* largebuf = MEMALLOC(sz  + strlen(tag1) + strlen(tag2) + strlen(tag3) + strlen(tag4),char);
718      largebuf[0] ='\0';      largebuf[0] ='\0';
719      char tmp[10];      char tmp[10];
720        __STRCAT(largebuf,tag1,tsz)
     __len = strlen(tag1);  
     __STRCAT(largebuf,tag1)  
     tsz +=  __len;  
   
721      for (i=numVTKNodesPerElement; i<=numGlobalCells*numVTKNodesPerElement; i+=numVTKNodesPerElement)      for (i=numVTKNodesPerElement; i<=numGlobalCells*numVTKNodesPerElement; i+=numVTKNodesPerElement)
722      {      {
723        sprintf(tmp,"%d\n", i);        sprintf(tmp,"%d\n", i);
724        __len=strlen(tmp);        __STRCAT(largebuf,tmp,tsz)
       __STRCAT(largebuf,tmp)  
       tsz+=__len;  
725      }      }
726        __STRCAT(largebuf,tag2,tsz)
727      __len=strlen(tag2);      largebuf[tsz] = '\0';
     __STRCAT(largebuf,tag2)  
     tsz+=__len;  
   
728      MPI_File_iwrite_shared(fh,largebuf, tsz,MPI_CHAR,&req);      MPI_File_iwrite_shared(fh,largebuf, tsz,MPI_CHAR,&req);
729      MPI_Wait(&req,&status);      MPI_Wait(&req,&status);
730        
731      // Re-using buffer!!      // re-using buffer!
732      largebuf[0] = '\0';      largebuf[0] = '\0';
733      tsz = 0;      tsz = 0;
734      __len = strlen(tag3);      __STRCAT(largebuf,tag3,tsz)
     __STRCAT(largebuf,tag3)  
     tsz+=__len;  
735      for (i=0; i<numGlobalCells; i++)      for (i=0; i<numGlobalCells; i++)
736      {      {
737        sprintf(tmp, "%d\n", cellType);        sprintf(tmp, "%d\n", cellType);
738        __len=strlen(tmp);        __STRCAT(largebuf,tmp,tsz)
739        __STRCAT(largebuf,tmp)      }
740        tsz+=__len;      __STRCAT(largebuf,tag4,tsz)
741      }      largebuf[tsz] = '\0';
     __len=strlen(tag4);  
     __STRCAT(largebuf,tag4)  
     tsz+=__len;  
   
742      MPI_File_iwrite_shared(fh,largebuf,tsz,MPI_CHAR,&req);      MPI_File_iwrite_shared(fh,largebuf,tsz,MPI_CHAR,&req);
743      MPI_Wait(&req,&status);      MPI_Wait(&req,&status);
744      MEMFREE(largebuf);      MEMFREE(largebuf);
# Line 772  void Finley_Mesh_saveVTK_MPIO(const char Line 746  void Finley_Mesh_saveVTK_MPIO(const char
746    
747    MPIO_DEBUG(" Done Writing Offsets & Types ")    MPIO_DEBUG(" Done Writing Offsets & Types ")
748    
749    // Write Point Data Header Tags    // Write Cell data header Tags
750    if( myRank == 0)    if(myRank == 0)
751    {    {
752      char header[600];      MPIO_DEBUG(" Writing Cell Data ...")
753      char tmpBuf[50];      if( write_celldata)
   
     if (write_pointdata)  
754      {      {
755        MPIO_DEBUG(" Writing Pointdata... ")        char tmpBuf[80];
756          char header[600];
757        // mark the active data arrays        // mark the active data arrays
758        bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;        bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
759        sprintf(header, "<PointData");        sprintf(tmpBuf, "<CellData");
760          strcat(header,tmpBuf);
761        for (i_data =0 ;i_data<num_data;++i_data)        for (i_data =0 ;i_data<num_data;++i_data)
762        {        {
763          if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])          if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data])
764          {          {
765            // if the rank == 0:   --> scalar data            // if the rank == 0:   --> scalar data
766            // if the rank == 1:   --> vector data            // if the rank == 1:   --> vector data
# Line 806  void Finley_Mesh_saveVTK_MPIO(const char Line 780  void Finley_Mesh_saveVTK_MPIO(const char
780              if (! set_vector)              if (! set_vector)
781              {              {
782                sprintf(tmpBuf," Vectors=\"%s\"",names_p[i_data]);                sprintf(tmpBuf," Vectors=\"%s\"",names_p[i_data]);
783                strcat(header,tmpBuf);            strcat(header,tmpBuf);
784                set_vector=TRUE;                set_vector=TRUE;
785              }              }
786              break;              break;
# Line 814  void Finley_Mesh_saveVTK_MPIO(const char Line 788  void Finley_Mesh_saveVTK_MPIO(const char
788              if (! set_tensor)              if (! set_tensor)
789              {              {
790                sprintf(tmpBuf," Tensors=\"%s\"",names_p[i_data]);                sprintf(tmpBuf," Tensors=\"%s\"",names_p[i_data]);
791                strcat(header,tmpBuf);            strcat(header,tmpBuf);
792                set_tensor=TRUE;                set_tensor=TRUE;
793              }              }
794              break;              break;
# Line 831  void Finley_Mesh_saveVTK_MPIO(const char Line 805  void Finley_Mesh_saveVTK_MPIO(const char
805      }      }
806    }    }
807    
808    // write actual data    // write actual data (collective)
809    if(write_pointdata)    if(write_celldata)
810    {    {
811      for (i_data =0 ;i_data<num_data;++i_data)      for (i_data =0 ;i_data<num_data;++i_data)
812      {      {
813        if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])        if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data])
814        {        {
815          numPointsPerSample = elements->ReferenceElement->numQuadNodes;          numPointsPerSample = elements->ReferenceElement->numQuadNodes;
816          rank = getDataPointRank(data_pp[i_data]);          rank = getDataPointRank(data_pp[i_data]);
# Line 871  void Finley_Mesh_saveVTK_MPIO(const char Line 845  void Finley_Mesh_saveVTK_MPIO(const char
845          if( myRank == 0)          if( myRank == 0)
846          {          {
847            char header[250];            char header[250];
848            sprintf(header, "<DataArray Name=\"%s\" type=\"Float32\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);            sprintf(header,"<DataArray Name=\"%s\" type=\"Float64\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);
849            MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);            MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);
850            MPI_Wait(&req,&status);            MPI_Wait(&req,&status);
851          }          }
         // write out the data  
         // if the number of required components is more than the number  
         // of actual components, pad with zeros  
852    
853          char tmpbuf[14];          // Write the actual data
854          char* largebuf = MEMALLOC(nCompReqd*numLocalNodes*14 + numLocal*nCompReqd + nCompReqd + 14,char);          char tmpbuf[15];
855            char* largebuf = MEMALLOC(nCompReqd*numLocalCells*15 + numLocalCells*nCompReqd + nCompReqd + 15,char);
856          largebuf[0] = '\0';          largebuf[0] = '\0';
         bool_t do_write=TRUE;  
857          size_t tsz = 0;          size_t tsz = 0;
858    
859          for(k=0;k < nodeCache.size;k++)          double sampleAvg[nComp];
860    
861            for (k=0; k<elementCache.size; k++)
862          {          {
863            i = nodeCache.values[k];            i = elementCache.values[k];
864    
865            if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)            values = getSampleData(data_pp[i_data], i);
866              // averaging over the number of points in the sample
867              for (n=0; n<nComp; n++)
868            {            {
869              if (mesh_p->Nodes->toReduced[i]>=0)              if (isExpanded(data_pp[i_data])) {
870              {                 rtmp = 0.;
871                switch(getFunctionSpaceType(data_pp[i_data]))                 for (j=0; j<numPointsPerSample; j++) rtmp += values[INDEX2(n,j,nComp)];
872                {                 sampleAvg[n] = rtmp/numPointsPerSample;
873                case FINLEY_DEGREES_OF_FREEDOM:              } else {
874                  values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);                 sampleAvg[n] = values[n];
                 break;  
               case FINLEY_REDUCED_DEGREES_OF_FREEDOM:  
                 values = getSampleData(data_pp[i_data],mesh_p->Nodes->reducedDegreeOfFreedom[i]);  
                 break;  
               case FINLEY_NODES:  
                 values = getSampleData(data_pp[i_data],i);  
                 break;  
               }  
               do_write=TRUE;  
             }  
             else  
             {  
               do_write=FALSE;  
875              }              }
876            }            }
877            else            // if the number of required components is more than the number
878              // of actual components, pad with zeros
879    
880              // probably only need to get shape of first element
881              // write the data different ways for scalar, vector and tensor
882              if (nCompReqd == 1)
883            {            {
884              do_write=TRUE;              sprintf(tmpbuf, " %e", sampleAvg[0]);
885              switch(getFunctionSpaceType(data_pp[i_data]))              __STRCAT(largebuf,tmpbuf,tsz)
             {  
             case FINLEY_DEGREES_OF_FREEDOM:  
               values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);  
               break;  
             case FINLEY_NODES:  
               values = getSampleData(data_pp[i_data],i);  
               break;  
             }  
886            }            }
887            // write the data different ways for scalar, vector and tensor            else if (nCompReqd == 3)
           if (do_write)  
888            {            {
889              if (nCompReqd == 1)              // write out the data
890                for (m=0; m<shape; m++)
891              {              {
892                sprintf(tmpbuf," %e", values[0]);                sprintf(tmpbuf, " %e", sampleAvg[m]);
893                __len=strlen(tmpbuf);                __STRCAT(largebuf,tmpbuf,tsz)
               __STRCAT(largebuf,tmpbuf)  
               tsz+=__len;  
894              }              }
895              else if (nCompReqd == 3)              for (m=0; m<nCompReqd-shape; m++)
896              {              {
897                for (m=0; m<shape; m++)                __STRCAT(largebuf,zero,tsz)
               {  
   
                 sprintf(tmpbuf," %e",values[m]);  
                 __len=strlen(tmpbuf);  
                 __STRCAT(largebuf,tmpbuf)  
                 tsz+=__len;  
               }  
               for (m=0; m<nCompReqd-shape; m++)  
               {  
                 __len=13;  
                 __STRCAT(largebuf,zero)  
                 tsz+=__len;  
               }  
898              }              }
899              else if (nCompReqd == 9)            }
900              else if (nCompReqd == 9)
901              {
902                // tensor data, so have a 3x3 matrix to output as a row
903                // of 9 data points
904                count = 0;
905                for (m=0; m<shape; m++)
906              {              {
907                // tensor data, so have a 3x3 matrix to output as a row                for (n=0; n<shape; n++)
               //  of 9 data points  
               count = 0;  
               for (m=0; m<shape; m++)  
908                {                {
909                  for (n=0; n<shape; n++)                  sprintf(tmpbuf, " %e", sampleAvg[count]);
910                  {                  __STRCAT(largebuf,tmpbuf,tsz)
911                    sprintf(tmpbuf," %e", values[count]);                  count++;
                   __len=strlen(tmpbuf);  
                   __STRCAT(largebuf,tmpbuf)  
                   tsz+=__len;  
                   count++;  
                 }  
                 for (n=0; n<3-shape; n++)  
                 {  
                   __len=13;  
                   __STRCAT(largebuf,zero)  
                   tsz+=__len;  
                 }  
912                }                }
913                for (m=0; m<3-shape; m++)                for (n=0; n<3-shape; n++)
914                {                {
915                  for (n=0; n<3; n++)                  __STRCAT(largebuf,zero,tsz)
                 {  
                   __len=13;  
                   __STRCAT(largebuf,zero)  
                   tsz+=__len;  
   
                 }  
916                }                }
917              }              }
918              __len=1;              for (m=0; m<3-shape; m++)
919              __STRCAT(largebuf,newline)                for (n=0; n<3; n++)
920              tsz+=1;                {
921                    __STRCAT(largebuf,zero,tsz)
922                  }
923            }            }
924              __STRCAT(largebuf,newline,tsz)
925          }          }
926          // Write out local data          largebuf[tsz] = '\0';
927          MPI_File_write_ordered(fh,largebuf,tsz,MPI_CHAR,&status);          MPI_File_write_ordered(fh,largebuf,tsz,MPI_CHAR,&status);
928          MEMFREE(largebuf);          MEMFREE(largebuf);
929          if( myRank == 0)          if( myRank == 0)
# Line 1000  void Finley_Mesh_saveVTK_MPIO(const char Line 932  void Finley_Mesh_saveVTK_MPIO(const char
932            MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);            MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);
933            MPI_Wait(&req,&status);            MPI_Wait(&req,&status);
934          }          }
935    
936        }        }
937      }      }
938      // Finish off with closing tag      // closing celldata tag
939      if(myRank == 0)      if(myRank == 0)
940      {      {
941        char* tag =  "</PointData>\n";        char* tag =  "</CellData>\n";
942        MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);        MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);
943        MPI_Wait(&req,&status);        MPI_Wait(&req,&status);
944      }      }
945      MPIO_DEBUG(" Done Writing Pointdata ")  
946        MPIO_DEBUG(" Done Writing Cell Data ")
947    }    }
   // end write_pointdata  
948    
949    // Write Cell data header Tags  
950    if(myRank == 0)    // Write Point Data Header Tags
951      if( myRank == 0)
952    {    {
953      if( write_celldata)      char header[600];
954        char tmpBuf[50];
955    
956        if (write_pointdata)
957      {      {
958        char tmpBuf[80];        MPIO_DEBUG(" Writing Pointdata... ")
       char header[600];  
959        // mark the active data arrays        // mark the active data arrays
960        bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;        bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
961        sprintf(tmpBuf, "<CellData");        sprintf(header, "<PointData");
       strcat(header,tmpBuf);  
962        for (i_data =0 ;i_data<num_data;++i_data)        for (i_data =0 ;i_data<num_data;++i_data)
963        {        {
964          if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data])          if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])
965          {          {
966            // if the rank == 0:   --> scalar data            // if the rank == 0:   --> scalar data
967            // if the rank == 1:   --> vector data            // if the rank == 1:   --> vector data
# Line 1046  void Finley_Mesh_saveVTK_MPIO(const char Line 981  void Finley_Mesh_saveVTK_MPIO(const char
981              if (! set_vector)              if (! set_vector)
982              {              {
983                sprintf(tmpBuf," Vectors=\"%s\"",names_p[i_data]);                sprintf(tmpBuf," Vectors=\"%s\"",names_p[i_data]);
984                  strcat(header,tmpBuf);
985                set_vector=TRUE;                set_vector=TRUE;
986              }              }
987              break;              break;
# Line 1053  void Finley_Mesh_saveVTK_MPIO(const char Line 989  void Finley_Mesh_saveVTK_MPIO(const char
989              if (! set_tensor)              if (! set_tensor)
990              {              {
991                sprintf(tmpBuf," Tensors=\"%s\"",names_p[i_data]);                sprintf(tmpBuf," Tensors=\"%s\"",names_p[i_data]);
992                  strcat(header,tmpBuf);
993                set_tensor=TRUE;                set_tensor=TRUE;
994              }              }
995              break;              break;
# Line 1069  void Finley_Mesh_saveVTK_MPIO(const char Line 1006  void Finley_Mesh_saveVTK_MPIO(const char
1006      }      }
1007    }    }
1008    
1009    // write actual data (collective)    // write actual data
1010    if(write_celldata)    if(write_pointdata)
1011    {    {
1012      for (i_data =0 ;i_data<num_data;++i_data)      for (i_data =0 ;i_data<num_data;++i_data)
1013      {      {
1014        if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data])        if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])
1015        {        {
1016          numPointsPerSample = elements->ReferenceElement->numQuadNodes;          numPointsPerSample = elements->ReferenceElement->numQuadNodes;
1017          rank = getDataPointRank(data_pp[i_data]);          rank = getDataPointRank(data_pp[i_data]);
# Line 1109  void Finley_Mesh_saveVTK_MPIO(const char Line 1046  void Finley_Mesh_saveVTK_MPIO(const char
1046          if( myRank == 0)          if( myRank == 0)
1047          {          {
1048            char header[250];            char header[250];
1049            sprintf(header,"<DataArray Name=\"%s\" type=\"Float32\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);            sprintf(header, "<DataArray Name=\"%s\" type=\"Float64\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);
1050            MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);            MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);
1051            MPI_Wait(&req,&status);            MPI_Wait(&req,&status);
1052          }          }
1053            // write out the data
1054            // if the number of required components is more than the number
1055            // of actual components, pad with zeros
1056    
1057          // Write the actual data */          char tmpbuf[15];
1058          char tmpbuf[14];          char* largebuf = MEMALLOC(nCompReqd*numLocalNodes*15 + numLocal*nCompReqd + nCompReqd + 15,char);
         char* largebuf = MEMALLOC(nCompReqd*numLocalCells*14 + numLocalCells*nCompReqd + nCompReqd + 14,char);  
1059          largebuf[0] = '\0';          largebuf[0] = '\0';
1060            bool_t do_write=TRUE;
1061          size_t tsz = 0;          size_t tsz = 0;
1062    
1063          double sampleAvg[nComp];          for(k=0;k < nodeCache.size;k++)
   
         for (k=0; i<elementCache.size; k++)  
1064          {          {
1065            i = elementCache.values[k];            i = nodeCache.values[k];
   
           values = getSampleData(data_pp[i_data], i);  
           // averaging over the number of points in the sample  
           for (k=0; k<nComp; k++)  
           {  
             rtmp = 0.;  
             for (j=0; j<numPointsPerSample; j++) rtmp += values[INDEX2(k,j,nComp)];  
             sampleAvg[k] = rtmp/numPointsPerSample;  
           }  
           // if the number of required components is more than the number  
           // of actual components, pad with zeros  
1066    
1067            // probably only need to get shape of first element            if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
           // write the data different ways for scalar, vector and tensor  
           if (nCompReqd == 1)  
           {  
             sprintf(tmpbuf, " %e", sampleAvg[0]);  
             __len=strlen(tmpbuf);  
             __STRCAT(largebuf,tmpbuf)  
             tsz+=__len;  
           }  
           else if (nCompReqd == 3)  
1068            {            {
1069              // write out the data              if (mesh_p->Nodes->toReduced[i]>=0)
             for (m=0; m<shape; m++)  
1070              {              {
1071                sprintf(tmpbuf, " %e", sampleAvg[m]);                switch(getFunctionSpaceType(data_pp[i_data]))
1072                __len=strlen(tmpbuf);                {
1073                __STRCAT(largebuf,tmpbuf)                case FINLEY_DEGREES_OF_FREEDOM:
1074                tsz+=__len;                  values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);
1075                    break;
1076                  case FINLEY_REDUCED_DEGREES_OF_FREEDOM:
1077                    values = getSampleData(data_pp[i_data],mesh_p->Nodes->reducedDegreeOfFreedom[i]);
1078                    break;
1079                  case FINLEY_NODES:
1080                    values = getSampleData(data_pp[i_data],i);
1081                    break;
1082                  }
1083                  do_write=TRUE;
1084              }              }
1085              for (m=0; m<nCompReqd-shape; m++)              else
1086              {              {
1087                tsz+=13;                do_write=FALSE;
               strcat(largebuf," 0.000000e+00");  
   
1088              }              }
1089            }            }
1090            else if (nCompReqd == 9)            else
1091            {            {
1092              // tensor data, so have a 3x3 matrix to output as a row              do_write=TRUE;
1093              // of 9 data points              switch(getFunctionSpaceType(data_pp[i_data]))
             count = 0;  
             for (m=0; m<shape; m++)  
1094              {              {
1095                for (n=0; n<shape; n++)              case FINLEY_DEGREES_OF_FREEDOM:
1096                  values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);
1097                  break;
1098                case FINLEY_NODES:
1099                  values = getSampleData(data_pp[i_data],i);
1100                  break;
1101                }
1102              }
1103              // write the data different ways for scalar, vector and tensor
1104              if (do_write)
1105              {
1106                if (nCompReqd == 1)
1107                {
1108                  sprintf(tmpbuf," %e", values[0]);
1109                  __STRCAT(largebuf,tmpbuf,tsz)
1110                }
1111                else if (nCompReqd == 3)
1112                {
1113                  for (m=0; m<shape; m++)
1114                {                {
                 sprintf(tmpbuf, " %e", sampleAvg[count]);  
                 __len=strlen(tmpbuf);  
                 __STRCAT(largebuf,tmpbuf)  
                 tsz+=__len;;  
   
1115    
1116                  count++;                  sprintf(tmpbuf," %e",values[m]);
1117                    __STRCAT(largebuf,tmpbuf,tsz)
1118                }                }
1119                for (n=0; n<3-shape; n++)                for (m=0; m<nCompReqd-shape; m++)
1120                {                {
1121                  __len=13;                  __STRCAT(largebuf,zero,tsz)
                 __STRCAT(largebuf,zero)  
                 tsz+=13;  
1122                }                }
1123              }              }
1124              for (m=0; m<3-shape; m++)              else if (nCompReqd == 9)
1125                for (n=0; n<3; n++)              {
1126                  // tensor data, so have a 3x3 matrix to output as a row
1127                  //  of 9 data points
1128                  count = 0;
1129                  for (m=0; m<shape; m++)
1130                  {
1131                    for (n=0; n<shape; n++)
1132                    {
1133                      sprintf(tmpbuf," %e", values[count]);
1134                      __STRCAT(largebuf,tmpbuf,tsz)
1135                      count++;
1136                    }
1137                    for (n=0; n<3-shape; n++)
1138                    {
1139                      __STRCAT(largebuf,zero,tsz)
1140                    }
1141                  }
1142                  for (m=0; m<3-shape; m++)
1143                {                {
1144                  __len=13;                  for (n=0; n<3; n++)
1145                  __STRCAT(largebuf,zero)                  {
1146                  tsz+=13;                    __STRCAT(largebuf,zero,tsz)
1147                    }
1148                }                }
1149                }
1150                __STRCAT(largebuf,newline,tsz)
1151            }            }
1152            __len=1;  
           __STRCAT(largebuf,newline)  
           tsz+=1;  
1153          }          }
1154            // Write out local data
1155    
1156            largebuf[tsz] = '\0';
1157          MPI_File_write_ordered(fh,largebuf,tsz,MPI_CHAR,&status);          MPI_File_write_ordered(fh,largebuf,tsz,MPI_CHAR,&status);
1158          MEMFREE(largebuf);          MEMFREE(largebuf);
1159          if( myRank == 0)          if( myRank == 0)
# Line 1207  void Finley_Mesh_saveVTK_MPIO(const char Line 1162  void Finley_Mesh_saveVTK_MPIO(const char
1162            MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);            MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);
1163            MPI_Wait(&req,&status);            MPI_Wait(&req,&status);
1164          }          }
   
1165        }        }
1166      }      }
1167      // closing celldata tag      // Finish off with closing tag
1168      if(myRank == 0)      if(myRank == 0)
1169      {      {
1170        char* tag =  "</CellData>\n";        char* tag =  "</PointData>\n";
1171        MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);        MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);
1172        MPI_Wait(&req,&status);        MPI_Wait(&req,&status);
1173      }      }
1174        MPIO_DEBUG(" Done Writing Pointdata ")
1175    }    }
1176      // end write_pointdata
1177    
1178    /* tag and bag... */    // tag and bag...  
1179    if (myRank == 0)    if (myRank == 0)
1180    {    {
1181      char *footer = "</Piece>\n</UnstructuredGrid>\n</VTKFile>";      char *footer = "</Piece>\n</UnstructuredGrid>\n</VTKFile>";
# Line 1234  void Finley_Mesh_saveVTK_MPIO(const char Line 1190  void Finley_Mesh_saveVTK_MPIO(const char
1190    MPI_Info_free(&infoHints);    MPI_Info_free(&infoHints);
1191  #undef MPIO_HINTS  #undef MPIO_HINTS
1192  #endif  #endif
   printf("\nTime: %f \n",  Paso_timer() - time0);  
1193    MPI_File_close(&fh);    MPI_File_close(&fh);
1194    MPIO_DEBUG(" ***** Exit saveVTK ***** ")    MPIO_DEBUG(" ***** Exit saveVTK ***** ")
1195  #undef __STRCAT  #undef __STRCAT
# Line 1248  void Finley_Mesh_saveVTK_MPIO(const char Line 1203  void Finley_Mesh_saveVTK_MPIO(const char
1203    
1204  void Finley_Mesh_saveVTK(const char * filename_p, Finley_Mesh *mesh_p, const dim_t num_data,char* *names_p, escriptDataC* *data_pp)  void Finley_Mesh_saveVTK(const char * filename_p, Finley_Mesh *mesh_p, const dim_t num_data,char* *names_p, escriptDataC* *data_pp)
1205  {  {
1206      #define NCOMP_MAX 9
1207    char error_msg[LenErrorMsg_MAX];    char error_msg[LenErrorMsg_MAX];
1208      double sampleAvg[NCOMP_MAX];
1209    /* if there is no mesh we just return */    /* if there is no mesh we just return */
   if (mesh_p==NULL) return;  
1210    
1211    int i, j, k, numVTKNodesPerElement,i_data,m, count, n, rank,shape, numPoints, cellType, numCells,    int i, j, k, numVTKNodesPerElement,i_data,m, count, n, rank,shape, numPoints, cellType, numCells,
1212    nDim, numPointsPerSample, nComp, nCompReqd;    nDim, numPointsPerSample, nComp, nCompReqd, NN;
   
1213    index_t j2;    index_t j2;
1214      int nodetype=FINLEY_DEGREES_OF_FREEDOM;
1215      int elementtype=FINLEY_UNKNOWN;
1216    double* values, rtmp;    double* values, rtmp;
1217    char elemTypeStr[32];    char elemTypeStr[32];
1218      FILE * fileHandle_p = NULL;
1219      bool_t do_write, *isCellCentered=NULL,write_celldata=FALSE,write_pointdata=FALSE;
1220      Finley_ElementFile* elements=NULL;
1221      ElementTypeId TypeId;
1222    
1223    printf("ddsafddfdafdf\n");
1224    /* open the file and check handle */    /* open the file and check handle */
1225      if (mesh_p==NULL) return;
1226    
1227    FILE * fileHandle_p = fopen(filename_p, "w");    fileHandle_p = fopen(filename_p, "w");
1228    if (fileHandle_p==NULL)    if (fileHandle_p==NULL)
1229    {    {
1230      sprintf(error_msg, "saveVTK: File %s could not be opened for writing.", filename_p);      sprintf(error_msg, "saveVTK: File %s could not be opened for writing.", filename_p);
# Line 1269  void Finley_Mesh_saveVTK(const char * fi Line 1232  void Finley_Mesh_saveVTK(const char * fi
1232      return;      return;
1233    }    }
1234    /* find the mesh type to be written */    /* find the mesh type to be written */
1235    int nodetype=FINLEY_DEGREES_OF_FREEDOM;    isCellCentered=TMPMEMALLOC(num_data,bool_t);
1236    int elementtype=FINLEY_UNKNOWN;  
1237    bool_t isCellCentered[num_data],write_celldata=FALSE,write_pointdata=FALSE;  
1238      if (Finley_checkPtr(isCellCentered)) {
1239         fclose(fileHandle_p);
1240         return;
1241      }
1242    for (i_data=0;i_data<num_data;++i_data)    for (i_data=0;i_data<num_data;++i_data)
1243    {    {
1244      if (! isEmpty(data_pp[i_data]))      if (! isEmpty(data_pp[i_data]))
# Line 1288  void Finley_Mesh_saveVTK(const char * fi Line 1255  void Finley_Mesh_saveVTK(const char * fi
1255          {          {
1256            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1257            fclose(fileHandle_p);            fclose(fileHandle_p);
1258              TMPMEMFREE(isCellCentered);
1259            return;            return;
1260          }          }
1261          isCellCentered[i_data]=FALSE;          isCellCentered[i_data]=FALSE;
# Line 1302  void Finley_Mesh_saveVTK(const char * fi Line 1270  void Finley_Mesh_saveVTK(const char * fi
1270          {          {
1271            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1272            fclose(fileHandle_p);            fclose(fileHandle_p);
1273              TMPMEMFREE(isCellCentered);
1274            return;            return;
1275          }          }
1276          isCellCentered[i_data]=FALSE;          isCellCentered[i_data]=FALSE;
# Line 1316  void Finley_Mesh_saveVTK(const char * fi Line 1285  void Finley_Mesh_saveVTK(const char * fi
1285          {          {
1286            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1287            fclose(fileHandle_p);            fclose(fileHandle_p);
1288              TMPMEMFREE(isCellCentered);
1289            return;            return;
1290          }          }
1291          isCellCentered[i_data]=FALSE;          isCellCentered[i_data]=FALSE;
# Line 1331  void Finley_Mesh_saveVTK(const char * fi Line 1301  void Finley_Mesh_saveVTK(const char * fi
1301            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1302            fclose(fileHandle_p);            fclose(fileHandle_p);
1303            return;            return;
1304              TMPMEMFREE(isCellCentered);
1305          }          }
1306          isCellCentered[i_data]=TRUE;          isCellCentered[i_data]=TRUE;
1307          break;          break;
# Line 1344  void Finley_Mesh_saveVTK(const char * fi Line 1315  void Finley_Mesh_saveVTK(const char * fi
1315          {          {
1316            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1317            fclose(fileHandle_p);            fclose(fileHandle_p);
1318              TMPMEMFREE(isCellCentered);
1319            return;            return;
1320          }          }
1321          isCellCentered[i_data]=TRUE;          isCellCentered[i_data]=TRUE;
# Line 1358  void Finley_Mesh_saveVTK(const char * fi Line 1330  void Finley_Mesh_saveVTK(const char * fi
1330          {          {
1331            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1332            fclose(fileHandle_p);            fclose(fileHandle_p);
1333              TMPMEMFREE(isCellCentered);
1334            return;            return;
1335          }          }
1336          isCellCentered[i_data]=TRUE;          isCellCentered[i_data]=TRUE;
# Line 1372  void Finley_Mesh_saveVTK(const char * fi Line 1345  void Finley_Mesh_saveVTK(const char * fi
1345          {          {
1346            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1347            fclose(fileHandle_p);            fclose(fileHandle_p);
1348              TMPMEMFREE(isCellCentered);
1349            return;            return;
1350          }          }
1351          isCellCentered[i_data]=TRUE;          isCellCentered[i_data]=TRUE;
# Line 1386  void Finley_Mesh_saveVTK(const char * fi Line 1360  void Finley_Mesh_saveVTK(const char * fi
1360          {          {
1361            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1362            fclose(fileHandle_p);            fclose(fileHandle_p);
1363              TMPMEMFREE(isCellCentered);
1364            return;            return;
1365          }          }
1366          isCellCentered[i_data]=TRUE;          isCellCentered[i_data]=TRUE;
# Line 1394  void Finley_Mesh_saveVTK(const char * fi Line 1369  void Finley_Mesh_saveVTK(const char * fi
1369          sprintf(error_msg,"saveVTK: Finley does not know anything about function space type %d",getFunctionSpaceType(data_pp[i_data]));          sprintf(error_msg,"saveVTK: Finley does not know anything about function space type %d",getFunctionSpaceType(data_pp[i_data]));
1370          Finley_setError(TYPE_ERROR,error_msg);          Finley_setError(TYPE_ERROR,error_msg);
1371          fclose(fileHandle_p);          fclose(fileHandle_p);
1372            TMPMEMFREE(isCellCentered);
1373          return;          return;
1374        }        }
1375        if (isCellCentered[i_data])        if (isCellCentered[i_data])
# Line 1417  void Finley_Mesh_saveVTK(const char * fi Line 1393  void Finley_Mesh_saveVTK(const char * fi
1393      numPoints = mesh_p->Nodes->numNodes;      numPoints = mesh_p->Nodes->numNodes;
1394    }    }
1395    if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;    if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;
   Finley_ElementFile* elements=NULL;  
1396    switch(elementtype)    switch(elementtype)
1397    {    {
1398    case FINLEY_ELEMENTS:    case FINLEY_ELEMENTS:
# Line 1437  void Finley_Mesh_saveVTK(const char * fi Line 1412  void Finley_Mesh_saveVTK(const char * fi
1412    {    {
1413      Finley_setError(SYSTEM_ERROR,"saveVTK: undefined element file");      Finley_setError(SYSTEM_ERROR,"saveVTK: undefined element file");
1414      fclose(fileHandle_p);      fclose(fileHandle_p);
1415        TMPMEMFREE(isCellCentered);
1416      return;      return;
1417    }    }
1418    /* map finley element type to VTK element type */    /* map finley element type to VTK element type */
1419    numCells = elements->numElements;    numCells = elements->numElements;
   ElementTypeId TypeId;  
1420    if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)    if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1421    {    {
1422      TypeId = elements->LinearReferenceElement->Type->TypeId;      TypeId = elements->LinearReferenceElement->Type->TypeId;
# Line 1550  void Finley_Mesh_saveVTK(const char * fi Line 1525  void Finley_Mesh_saveVTK(const char * fi
1525      sprintf(error_msg, "saveVTK: Element type %s is not supported by VTK",elements->ReferenceElement->Type->Name);      sprintf(error_msg, "saveVTK: Element type %s is not supported by VTK",elements->ReferenceElement->Type->Name);
1526      Finley_setError(VALUE_ERROR,error_msg);      Finley_setError(VALUE_ERROR,error_msg);
1527      fclose(fileHandle_p);      fclose(fileHandle_p);
1528        TMPMEMFREE(isCellCentered);
1529      return;      return;
1530    }    }
1531    /* xml header */    /* xml header */
# Line 1572  void Finley_Mesh_saveVTK(const char * fi Line 1548  void Finley_Mesh_saveVTK(const char * fi
1548    * the reason for this if statement is explained in the long comment below    * the reason for this if statement is explained in the long comment below
1549    */    */
1550    nDim = mesh_p->Nodes->numDim;    nDim = mesh_p->Nodes->numDim;
1551    fprintf(fileHandle_p, "<DataArray NumberOfComponents=\"%d\" type=\"Float32\" format=\"ascii\">\n",MAX(3,nDim));    fprintf(fileHandle_p, "<DataArray NumberOfComponents=\"%d\" type=\"Float64\" format=\"ascii\">\n",MAX(3,nDim));
1552    /* vtk/mayavi doesn't like 2D data, it likes 3D data with a degenerate    /* vtk/mayavi doesn't like 2D data, it likes 3D data with a degenerate
1553    * third dimension to handle 2D data (like a sheet of paper).  So, if    * third dimension to handle 2D data (like a sheet of paper).  So, if
1554    * nDim is 2, we have to append zeros to the array to get this third    * nDim is 2, we have to append zeros to the array to get this third
# Line 1608  void Finley_Mesh_saveVTK(const char * fi Line 1584  void Finley_Mesh_saveVTK(const char * fi
1584    
1585    /* write out the DataArray element for the connectivity */    /* write out the DataArray element for the connectivity */
1586    
1587    int NN = elements->ReferenceElement->Type->numNodes;    NN = elements->ReferenceElement->Type->numNodes;
1588    fprintf(fileHandle_p, "<Cells>\n");    fprintf(fileHandle_p, "<Cells>\n");
1589    fprintf(fileHandle_p, "<DataArray Name=\"connectivity\" type=\"Int32\" format=\"ascii\">\n");    fprintf(fileHandle_p, "<DataArray Name=\"connectivity\" type=\"Int32\" format=\"ascii\">\n");
1590    
# Line 1720  void Finley_Mesh_saveVTK(const char * fi Line 1696  void Finley_Mesh_saveVTK(const char * fi
1696            sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);            sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);
1697            Finley_setError(VALUE_ERROR,error_msg);            Finley_setError(VALUE_ERROR,error_msg);
1698            fclose(fileHandle_p);            fclose(fileHandle_p);
1699              TMPMEMFREE(isCellCentered);
1700            return;            return;
1701          }          }
1702        }        }
# Line 1746  void Finley_Mesh_saveVTK(const char * fi Line 1723  void Finley_Mesh_saveVTK(const char * fi
1723            {            {
1724              Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");              Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");
1725              fclose(fileHandle_p);              fclose(fileHandle_p);
1726                TMPMEMFREE(isCellCentered);
1727              return;              return;
1728            }            }
1729            nCompReqd = 3;            nCompReqd = 3;
# Line 1757  void Finley_Mesh_saveVTK(const char * fi Line 1735  void Finley_Mesh_saveVTK(const char * fi
1735            {            {
1736              Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");              Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
1737              fclose(fileHandle_p);              fclose(fileHandle_p);
1738                TMPMEMFREE(isCellCentered);
1739              return;              return;
1740            }            }
1741            nCompReqd = 9;            nCompReqd = 9;
1742          }          }
1743          fprintf(fileHandle_p, "<DataArray Name=\"%s\" type=\"Float32\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);          fprintf(fileHandle_p, "<DataArray Name=\"%s\" type=\"Float64\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);
1744    
         double sampleAvg[nComp];  
1745          for (i=0; i<numCells; i++)          for (i=0; i<numCells; i++)
1746          {          {
1747            values = getSampleData(data_pp[i_data], i);            values = getSampleData(data_pp[i_data], i);
1748            /* averaging over the number of points in the sample */            /* averaging over the number of points in the sample */
1749            for (k=0; k<nComp; k++)            for (k=0; k<MIN(nComp,NCOMP_MAX); k++)
1750            {            {
1751              rtmp = 0.;              if (isExpanded(data_pp[i_data])) {
1752              for (j=0; j<numPointsPerSample; j++) rtmp += values[INDEX2(k,j,nComp)];                 rtmp = 0.;
1753              sampleAvg[k] = rtmp/numPointsPerSample;                 for (j=0; j<numPointsPerSample; j++) rtmp += values[INDEX2(k,j,nComp)];
1754                   sampleAvg[k] = rtmp/numPointsPerSample;
1755                } else {
1756                   sampleAvg[k] = values[k];
1757                }
1758    
1759            }            }
1760            /* if the number of required components is more than the number            /* if the number of required components is more than the number
1761            * of actual components, pad with zeros            * of actual components, pad with zeros
# Line 1854  void Finley_Mesh_saveVTK(const char * fi Line 1837  void Finley_Mesh_saveVTK(const char * fi
1837            sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);            sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);
1838            Finley_setError(VALUE_ERROR,error_msg);            Finley_setError(VALUE_ERROR,error_msg);
1839            fclose(fileHandle_p);            fclose(fileHandle_p);
1840              TMPMEMFREE(isCellCentered);
1841            return;            return;
1842          }          }
1843        }        }
# Line 1880  void Finley_Mesh_saveVTK(const char * fi Line 1864  void Finley_Mesh_saveVTK(const char * fi
1864            {            {
1865              Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");              Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");
1866              fclose(fileHandle_p);              fclose(fileHandle_p);
1867                TMPMEMFREE(isCellCentered);
1868              return;              return;
1869            }            }
1870            nCompReqd = 3;            nCompReqd = 3;
# Line 1891  void Finley_Mesh_saveVTK(const char * fi Line 1876  void Finley_Mesh_saveVTK(const char * fi
1876            {            {
1877              Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");              Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
1878              fclose(fileHandle_p);              fclose(fileHandle_p);
1879                TMPMEMFREE(isCellCentered);
1880              return;              return;
1881            }            }
1882            nCompReqd = 9;            nCompReqd = 9;
1883          }          }
1884          fprintf(fileHandle_p, "<DataArray Name=\"%s\" type=\"Float32\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);          fprintf(fileHandle_p, "<DataArray Name=\"%s\" type=\"Float64\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);
1885          /* write out the data */          /* write out the data */
1886          /* if the number of required components is more than the number          /* if the number of required components is more than the number
1887          * of actual components, pad with zeros          * of actual components, pad with zeros
1888          */          */
1889          bool_t do_write=TRUE;          do_write=TRUE;
1890          for (i=0; i<mesh_p->Nodes->numNodes; i++)          for (i=0; i<mesh_p->Nodes->numNodes; i++)
1891          {          {
1892            if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)            if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
# Line 1984  void Finley_Mesh_saveVTK(const char * fi Line 1970  void Finley_Mesh_saveVTK(const char * fi
1970    fprintf(fileHandle_p, "</VTKFile>\n");    fprintf(fileHandle_p, "</VTKFile>\n");
1971    /* close the file */    /* close the file */
1972    fclose(fileHandle_p);    fclose(fileHandle_p);
1973      TMPMEMFREE(isCellCentered);
1974    return;    return;
1975  }  }
1976  #endif  #endif

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

  ViewVC Help
Powered by ViewVC 1.1.26