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

Diff of /trunk-mpi-branch/finley/src/Mesh_saveVTK.c

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

revision 1222 by ksteube, Tue May 15 03:23:17 2007 UTC revision 1223 by gross, Fri Aug 3 02:40:39 2007 UTC
# Line 14  Line 14 
14  /**************************************************************/  /**************************************************************/
15    
16  /*   writes data and mesh in a vtk file */  /*   writes data and mesh in a vtk file */
17    /*   nodal data needs to be given on FINLEY_NODES or FINLEY_REDUCED_NODES */
18    
19  /**************************************************************/  /**************************************************************/
20    
21  /*   Author: Paul Cochrane, cochrane@esscc.uq.edu.au */  /*   Author: Paul Cochrane, cochrane@esscc.uq.edu.au */
22  /*   MPI-IO version: Derick Hawcroft, d.hawcroft@uq.edu.au         */  /*   MPI version: Derick Hawcroft, d.hawcroft@uq.edu.au  */
23    
24  /*   Version: $Id$ */  /*   Version: $Id$ */
25    
# Line 29  Line 30 
30  #include "Assemble.h"  #include "Assemble.h"
31  #include "vtkCellType.h"  /* copied from vtk source directory !!! */  #include "vtkCellType.h"  /* copied from vtk source directory !!! */
32    
33  /*  #define LEN_PRINTED_INT_FORMAT (9+1)
34   MPI version notes:  #define INT_FORMAT "%d "
35    #define INT_NEWLINE_FORMAT "%d\n"
36   ******************************************************************************  #define FLOAT_SCALAR_FORMAT "%15.10e\n"
37   ***                                                                       ****  #define FLOAT_VECTOR_FORMAT "%15.10e %15.10e %15.10e\n"
38   *** WARNING: Won't work for meshes with peridodic boundary conditions yet ****  #define FLOAT_TENSOR_FORMAT "%15.10e %15.10e %15.10e %15.10e %15.10e %15.10e %15.10e %15.10e %15.10e\n"
39   ***                                                                       ****    #define LEN_PRINTED_FLOAT_SCALAR_FORMAT (15+1)
40   ******************************************************************************  #define LEN_PRINTED_FLOAT_VECTOR_FORMAT (3*(15+1)+1)
41    #define LEN_PRINTED_FLOAT_TENSOR_FORMAT (9*(15+1)+1)
42   In this version, the rank==0 process writes *all* opening and closing  #define NEWLINE "\n"
43   XML tags.  #define LEN_TMP_BUFFER LEN_PRINTED_FLOAT_TENSOR_FORMAT+(MAX_numNodes*LEN_PRINTED_INT_FORMAT+1)+2
44   Individual process data is copied to a buffer before being written  #define NCOMP_MAX 9
45   out. The  routines are collectively called and will be called in the natural  #define __STRCAT(dest,chunk,dest_in_use)  \
46   ordering i.e 0 to maxProcs-1.  {                  \
47      strcpy(&dest[dest_in_use], chunk); \
48   Notable Notables:    dest_in_use+=strlen(chunk); \
  the struct localIndexCache: stores local domain indices for faster  reference  
 */  
   
 #ifdef PASO_MPI  
   
   
 //#define MPIO_HINTS  
   
   
   
 #define MPIO_DEBUG(str) \  
 { \  
     if(myRank == 0) \  
     printf("==== MPI-IO => %s \n", str); \  
49  }  }
50    
51  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(const char * filename_p,
52                             Finley_Mesh *mesh_p,
53                             const dim_t num_data,
54                             char* *names_p,
55                             escriptDataC* *data_pp)
56  {  {
57    int    numPoints,    char error_msg[LenErrorMsg_MAX], *txt_buffer=NULL, tmp_buffer[LEN_TMP_BUFFER];
58    numCells = -1,    double sampleAvg[NCOMP_MAX], *values, rtmp;
59               myRank,comm,gsize,    size_t len_txt_buffer, max_len_names, txt_buffer_in_use;
60               numLocal,    FILE * fileHandle_p = NULL;
61               nDim,    int mpi_size, i_data, i,j , cellType;
62               shape;    dim_t nDim, globalNumPoints, numCells, globalNumCells, numVTKNodesPerElement, myNumPoints, numPointsPerSample, rank, nComp, nCompReqd, shape, NN, numCellFactor, myNumCells;
63    size_t __n;    bool_t do_write, *isCellCentered=NULL,write_celldata=FALSE,write_pointdata=FALSE;
64    int i,j,k,m,n,count;    bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
65    int numGlobalCells = 0;    index_t myFirstNode, myLastNode, *globalNodeIndex, k, *node_index, myFirstCell;
66    index_t  *nodesGlobal=NULL;   // used to get the connectivity  right for VTK    #ifdef PASO_MPI
67      int ierr;
   /* variables associatted with write_celldata/pointdata */  
   int numPointsPerSample,  
   nComp,  
   nCompReqd;  
   double* values, rtmp;  
   
   // Local element info (for debugging)  
   size_t numLocalCells,  
   numInternalCells,  
   numBoundaryCells;  
   
   int rank;  
   
68    int amode = MPI_MODE_CREATE | MPI_MODE_WRONLY |  MPI_MODE_SEQUENTIAL;    int amode = MPI_MODE_CREATE | MPI_MODE_WRONLY |  MPI_MODE_SEQUENTIAL;
69      MPI_File mpi_fileHandle_p;
70    comm   = mesh_p->Nodes->MPIInfo->comm;    MPI_Status mpi_status;
71    myRank = mesh_p->Nodes->MPIInfo->rank;    MPI_Request mpi_req;
72    gsize  = mesh_p->Nodes->MPIInfo->size;    MPI_Info mpi_info=MPI_INFO_NULL;
73      #endif
74    MPI_File fh;    Paso_MPI_rank my_mpi_rank;
75    MPI_Status status;    int nodetype=FINLEY_NODES;
   MPI_Request req;  
   MPI_Info infoHints;  
   
   char error_msg[LenErrorMsg_MAX];  
   
   int i_data;  
   
   int nodetype=FINLEY_DEGREES_OF_FREEDOM;  
76    int elementtype=FINLEY_UNKNOWN;    int elementtype=FINLEY_UNKNOWN;
   bool_t isCellCentered[num_data],write_celldata=FALSE,write_pointdata=FALSE;  
   
   ElementTypeId TypeId;  
   
   int numVTKNodesPerElement;  
   int cellType;  
77    char elemTypeStr[32];    char elemTypeStr[32];
78      Finley_NodeMapping *nodeMapping=NULL;
79    Finley_ElementFile* elements=NULL;    Finley_ElementFile* elements=NULL;
80      ElementTypeId TypeId;
   
   // Local node info  
   int numInternalNodes,  
   numLocalNodes,  
   numBoundaryNodes,  
   localDOF;  // equals to  (DOF of Internal Nodes) +  (DOF of Boundary Nodes) of local domain  
   
   
   nDim  = mesh_p->Nodes->numDim;  
   
 #ifdef MPIO_HINTS  
   MPI_Info_create(&infoHints);  
   //  MPI_Info_set(infoHints, "striping_unit",        "424288");  
   //  MPI_Info_set(infoHints, "striping_factor",      "16");  
   //  MPI_Info_set(infoHints, "collective_buffering", "true");  
   //  MPI_Info_set(infoHints, "cb_block_size",        "131072");  
   //  MPI_Info_set(infoHints, "cb_buffer_size",       "1048567");  
   //  MPI_Info_set(infoHints, "cb_nodes",             "8");  
   //    MPI_Info_set(infoHints, "access_style", "write_once, sequential");  
   
   //XFS only  
   //   MPI_Info_set(infoHints, "direct_write",          "true");  
 #else  
   infoHints = MPI_INFO_NULL;  
 #endif  
   
   // Holds a local node/element index into the global array  
   struct localIndexCache  
   {  
     index_t *values;  
     int size;  
   };  
   
   struct localIndexCache nodeCache,  
           elementCache;  
   
   // Collective Call  
   MPI_File_open(mesh_p->Nodes->MPIInfo->comm, (char*)filename_p, amode,infoHints, &fh);  
   MPI_File_set_view(fh,MPI_DISPLACEMENT_CURRENT,MPI_CHAR, MPI_CHAR, "native" , infoHints);  
   
   MPIO_DEBUG(" ***** Enter saveVTK ******")  
   
   for (i_data=0;i_data<num_data;++i_data)  
   {  
     if (! isEmpty(data_pp[i_data]) )  
     {  
       switch(getFunctionSpaceType(data_pp[i_data]))  
       {  
       case FINLEY_DEGREES_OF_FREEDOM:  
         nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;  
         if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)  
         {  
           elementtype=FINLEY_ELEMENTS;  
         }  
         else  
         {  
           Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");  
           return ;  
         }  
         isCellCentered[i_data]=FALSE;  
         break;  
       case FINLEY_REDUCED_DEGREES_OF_FREEDOM:  
         nodetype = FINLEY_REDUCED_DEGREES_OF_FREEDOM;  
         if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)  
         {  
           elementtype=FINLEY_ELEMENTS;  
         }  
         else  
         {  
           Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");  
           return;  
         }  
         isCellCentered[i_data]=FALSE;  
         break;  
       case FINLEY_NODES:  
         nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;  
         if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)  
         {  
           elementtype=FINLEY_ELEMENTS;  
         }  
         else  
         {  
           Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");  
           return;  
         }  
         isCellCentered[i_data]=FALSE;  
         break;  
       case FINLEY_ELEMENTS:  
       case FINLEY_REDUCED_ELEMENTS:  
         nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;  
         if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)  
         {  
           elementtype=FINLEY_ELEMENTS;  
         }  
         else  
         {  
           Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");  
           return;  
         }  
         isCellCentered[i_data]=TRUE;  
         break;  
       case FINLEY_FACE_ELEMENTS:  
       case FINLEY_REDUCED_FACE_ELEMENTS:  
         nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;  
         if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_FACE_ELEMENTS)  
         {  
           elementtype=FINLEY_FACE_ELEMENTS;  
         }  
         else  
         {  
           Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");  
           return;  
   
         }  
         isCellCentered[i_data]=TRUE;  
         break;  
       case FINLEY_POINTS:  
         nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;  
         if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_POINTS)  
         {  
           elementtype=FINLEY_POINTS;  
         }  
         else  
         {  
           Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");  
           return;  
   
         }  
         isCellCentered[i_data]=TRUE;  
         break;  
       case FINLEY_CONTACT_ELEMENTS_1:  
       case FINLEY_REDUCED_CONTACT_ELEMENTS_1:  
         nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;  
         if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1)  
         {  
           elementtype=FINLEY_CONTACT_ELEMENTS_1;  
         }  
         else  
         {  
           Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");  
           return;  
   
         }  
         isCellCentered[i_data]=TRUE;  
         break;  
       case FINLEY_CONTACT_ELEMENTS_2:  
       case FINLEY_REDUCED_CONTACT_ELEMENTS_2:  
         nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;  
         if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1)  
         {  
           elementtype=FINLEY_CONTACT_ELEMENTS_1;  
         }  
         else  
         {  
           Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");  
           return;  
   
         }  
         isCellCentered[i_data]=TRUE;  
         break;  
       default:  
         sprintf(error_msg,"saveVTK: Finley does not know anything about function space type %d",getFunctionSpaceType(data_pp[i_data]));  
         Finley_setError(TYPE_ERROR,error_msg);  
         return;  
   
       }  
   
       if (isCellCentered[i_data])  
       {  
         write_celldata=TRUE;  
       }  
       else  
       {  
         write_pointdata=TRUE;  
       }  
     }  
   }  
   
   Finley_NodeDistribution *dist;  
   if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)  
   {  
     dist = mesh_p->Nodes->reducedDegreeOfFreedomDistribution;  
   }  
   else  
   {  
     dist = mesh_p->Nodes->degreeOfFreedomDistribution;  
   }  
   
   numInternalNodes = dist->numInternal;  
   numBoundaryNodes = dist->numBoundary;  
   
   localDOF =  dist->numLocal;  
   
   numPoints        = dist->numGlobal;  
   
   if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;  
   switch(elementtype)  
   {  
   case FINLEY_ELEMENTS:  
     elements=mesh_p->Elements;  
     break;  
   case FINLEY_FACE_ELEMENTS:  
     elements=mesh_p->FaceElements;  
     break;  
   case FINLEY_POINTS:  
     elements=mesh_p->Points;  
     break;  
   case FINLEY_CONTACT_ELEMENTS_1:  
     elements=mesh_p->ContactElements;  
     break;  
   }  
   if (elements==NULL)  
   {  
     Finley_setError(SYSTEM_ERROR,"saveVTK: undefined element file");  
     return ;  
   }  
   
   numCells =  elements->numElements;  
   numGlobalCells = elements->elementDistribution->vtxdist[gsize];  
   numLocalCells    = elements->elementDistribution->numLocal;  
   numInternalCells = elements->elementDistribution->numInternal;  
   numBoundaryCells = elements->elementDistribution->numBoundary;  
   
   if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM || nodetype==FINLEY_REDUCED_NODES)  
   {  
     TypeId = elements->LinearReferenceElement->Type->TypeId;  
   }  
   else  
   {  
     TypeId = elements->ReferenceElement->Type->TypeId;  
   }  
   
   switch(TypeId)  
   {  
   case Point1:  
   case Line2Face:  
   case Line3Face:  
   case Point1_Contact:  
   case Line2Face_Contact:  
   case Line3Face_Contact:  
     cellType = VTK_VERTEX;  
     numVTKNodesPerElement = 1;  
     strcpy(elemTypeStr, "VTK_VERTEX");  
     break;  
   
   case Line2:  
   case Tri3Face:  
   case Rec4Face:  
   case Line2_Contact:  
   case Tri3_Contact:  
   case Tri3Face_Contact:  
   case Rec4Face_Contact:  
     cellType = VTK_LINE;  
     numVTKNodesPerElement = 2;  
     strcpy(elemTypeStr, "VTK_LINE");  
     break;  
   
   case Tri3:  
   case Tet4Face:  
   case Tet4Face_Contact:  
     cellType = VTK_TRIANGLE;  
     numVTKNodesPerElement = 3;  
     strcpy(elemTypeStr, "VTK_TRIANGLE");  
     break;  
   
   case Rec4:  
   case Hex8Face:  
   case Rec4_Contact:  
   case Hex8Face_Contact:  
     cellType = VTK_QUAD;  
     numVTKNodesPerElement = 4;  
     strcpy(elemTypeStr, "VTK_QUAD");  
     break;  
   
   case Tet4:  
     cellType = VTK_TETRA;  
     numVTKNodesPerElement = 4;  
     strcpy(elemTypeStr, "VTK_TETRA");  
     break;  
   
   case Hex8:  
     cellType = VTK_HEXAHEDRON;  
     numVTKNodesPerElement = 8;  
     strcpy(elemTypeStr, "VTK_HEXAHEDRON");  
     break;  
   
   case Line3:  
   case Tri6Face:  
   case Rec8Face:  
   case Line3_Contact:  
   case Tri6Face_Contact:  
   case Rec8Face_Contact:  
     cellType = VTK_QUADRATIC_EDGE;  
     numVTKNodesPerElement = 3;  
     strcpy(elemTypeStr, "VTK_QUADRATIC_EDGE");  
     break;  
   
   case Tri6:  
   case Tet10Face:  
   case Tri6_Contact:  
   case Tet10Face_Contact:  
     cellType = VTK_QUADRATIC_TRIANGLE;  
     numVTKNodesPerElement = 6;  
     strcpy(elemTypeStr, "VTK_QUADRATIC_TRIANGLE");  
     break;  
   
   case Rec8:  
   case Hex20Face:  
   case Rec8_Contact:  
   case Hex20Face_Contact:  
     cellType = VTK_QUADRATIC_QUAD;  
     numVTKNodesPerElement = 8;  
     strcpy(elemTypeStr, "VTK_QUADRATIC_QUAD");  
     break;  
   
   case Tet10:  
     cellType = VTK_QUADRATIC_TETRA;  
     numVTKNodesPerElement = 10;  
     strcpy(elemTypeStr, "VTK_QUADRATIC_TETRA");  
     break;  
   
   case Hex20:  
     cellType = VTK_QUADRATIC_HEXAHEDRON;  
     numVTKNodesPerElement = 20;  
     strcpy(elemTypeStr, "VTK_QUADRATIC_HEXAHEDRON");  
     break;  
   
   default:  
     sprintf(error_msg, "saveVTK: Element type %s is not supported by VTK",elements->ReferenceElement->Type->Name);  
     Finley_setError(VALUE_ERROR,error_msg);  
     return;  
   }  
   
   /* Write XML Header */  
   if(myRank == 0)  
   {  
     char header[400];  
   
     sprintf(header,"<?xml version=\"1.0\"?>\n" \  
             "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\">\n" \  
             "<UnstructuredGrid>\n" \  
             "<Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n" \  
             "<Points>\n" \  
             "<DataArray NumberOfComponents=\"%d\" type=\"Float64\" format=\"ascii\">\n"  
             ,numPoints,numGlobalCells,MAX(3,nDim));  
   
   
     MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);  
     MPI_Wait(&req,&status);  
   }  
   
   MPIO_DEBUG(" Writing Coordinate Points... ")  
   
   numLocalNodes=localDOF;  
   
   //  we will be writing values which vary from 13-15 chars hence the strlen()  
   char* largebuf = MEMALLOC( numLocalNodes*15*nDim + numLocalNodes*2 + 1 ,char);  
   largebuf[0] = '\0';  
   char tmpbuf[15];  
   int tsz=0;  
   int numNodesOutput=0;  
   index_t pos=0;  
   
   index_t *vtxdist = NULL, *DOFNodes=NULL,*forwardBuffer=NULL,*backwardBuffer=NULL;  
   
   DOFNodes   = MEMALLOC(mesh_p->Nodes->numNodes,index_t);  
   nodeCache.values = MEMALLOC( numLocalNodes, index_t);  
   index_t bc_pos = 0;  
   
   // Custom string concat:  avoids expensive strlen(3) call by strcat(3)  
   // Note the implicit assumption on the variable "tsz"  
   int __len,__j;  
   char  *zero = "0.000000e+00 ";  
   char  *newline = "\n";  
     
 #define __STRCAT(dest,chunk,tsz)  \  
 {                  \  
    __len = strlen(chunk); \  
    __j = -1;      \  
    while(__j++ < __len)  \  
     *(dest+tsz+__j)=*(chunk+__j); \  
    tsz+=__len;              \  
 }  
81        
82    // Loop over all nodes      
83    for (i = 0; i < mesh_p->Nodes->numNodes; i++)    /****************************************/
84    {    /*                                      */
85      // This is the bit that will break for periodic BCs because it assumes that there is a one to one    /*       tags in the vtk file           */
86      // correspondance between nodes and Degrees of freedom  
87      //TODO: handle periodic BC's    char* tags_header="<?xml version=\"1.0\"?>\n" \
88      DOFNodes[mesh_p->Nodes->degreeOfFreedom[i]] = i;                      "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\">\n" \
89                        "<UnstructuredGrid>\n" \
90      // Is this node local to the domain ?                      "<Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n" \
91      if( mesh_p->Nodes->degreeOfFreedom[i] < localDOF )                      "<Points>\n" \
92      {                      "<DataArray NumberOfComponents=\"%d\" type=\"Float64\" format=\"ascii\">\n";
93        for (j = 0; j < nDim; j++)    char *tag_End_DataArray = "</DataArray>\n";
94        {    char* tag_End_PointData =  "</PointData>\n";
95          sprintf(tmpbuf,"%e ", mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)] );    char* tag_End_CellData =  "</CellData>\n";
96          __STRCAT(largebuf,tmpbuf,tsz)    char *footer = "</Piece>\n</UnstructuredGrid>\n</VTKFile>";
97        }    char* tags_End_Points_and_Start_Conn = "</DataArray>\n</Points>\n<Cells>\n<DataArray Name=\"connectivity\" type=\"Int32\" format=\"ascii\">\n" ;
98        for (k=0; k<3-nDim; k++)    char* tags_End_Conn_and_Start_Offset = "</DataArray>\n<DataArray Name=\"offsets\" type=\"Int32\" format=\"ascii\">\n";
99        {    char* tags_End_Offset_and_Start_Type = "</DataArray>\n<DataArray Name=\"types\" type=\"UInt8\" format=\"ascii\">\n";
100          __STRCAT(largebuf,zero,tsz)    char* tag_Float_DataArray="<DataArray Name=\"%s\" type=\"Float64\" NumberOfComponents=\"%d\" format=\"ascii\">\n";
101        }    char *tags_End_Type_And_Cells = "</DataArray>\n</Cells>\n";
       __STRCAT(largebuf,newline,tsz)  
       nodeCache.values[numNodesOutput++]=i;  
     }  
   }  
   
   nodeCache.size=numNodesOutput;  
102    
103    largebuf[tsz] = '\0';    int VTK_QUADRATIC_HEXAHEDRON_INDEX[] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 16, 17, 18, 19, 12, 13, 14, 15 };
104    MPI_File_write_ordered(fh, largebuf,tsz, MPI_CHAR, &status);    /* if there is no mesh we just return */
105    MEMFREE(largebuf);    if (mesh_p==NULL) return;
   
   nodesGlobal = MEMALLOC(mesh_p->Nodes->numNodes ,index_t);  
   
   // form distribution info on who output which nodes  
   vtxdist = MEMALLOC( gsize+1, index_t );  
   vtxdist[0]=0;  
   MPI_Allgather(&numNodesOutput,1,MPI_INT,vtxdist+1,1,MPI_INT,comm);  
   for( i=0; i<gsize; i++ )  
     vtxdist[i+1]+=vtxdist[i];  
   
   // will not work for periodic boundary conditions  
   // calculate the local nodes file positions  
   pos = 0;  
   for( i=0; i<mesh_p->Nodes->numNodes; i++ )  
   {  
     if( mesh_p->Nodes->degreeOfFreedom[i]< localDOF )  
     {  
       nodesGlobal[i] = vtxdist[myRank] + pos++;  
     }  
     else  
       nodesGlobal[i] = -1;  
   }  
106    
107    // communicate the local Nodes file position to the interested parties    my_mpi_rank = mesh_p->Nodes->MPIInfo->rank;
108    // send local info    mpi_size  = mesh_p->Nodes->MPIInfo->size;
109    forwardBuffer = MEMALLOC( mesh_p->Nodes->numNodes, index_t );    nDim = mesh_p->Nodes->numDim;
   for( n=0; n < dist->numNeighbours; n++ )  
   {  
     if(  dist->edges[n]->numForward)  
     {  
       for( i=0; i < dist->edges[n]->numForward; i++ )  
         forwardBuffer[i] = nodesGlobal[DOFNodes[dist->edges[n]->indexForward[i] ]];  
       Paso_CommBuffer_pack( mesh_p->Nodes->CommBuffer, dist->neighbours[n], NULL, forwardBuffer, sizeof(index_t), 0 );  
       Paso_CommBuffer_send( mesh_p->Nodes->CommBuffer, dist->neighbours[n], sizeof(index_t) );  
     }  
   }  
   // receive external info  
   backwardBuffer = MEMALLOC( mesh_p->Nodes->numNodes, index_t );  
   for( n=0; n < dist->numNeighbours; n++ )  
   {  
     if( dist->edges[n]->numBackward )  
     {  
       Paso_CommBuffer_recv(mesh_p->Nodes->CommBuffer, dist->neighbours[n], sizeof(index_t));  
       Paso_CommBuffer_unpack(mesh_p->Nodes->CommBuffer, dist->neighbours[n], NULL, backwardBuffer, sizeof(index_t), 0 );  
       /* TODO: voodoo to handle perdiodic  BC's */  
       for( i=0; i<dist->edges[n]->numBackward; i++ )  
         nodesGlobal[DOFNodes[dist->edges[n]->indexBackward[i] ]] = backwardBuffer[i];  
     }  
   }  
     
110    
111        if (! (nDim ==2) || (nDim == 3) ) {
112    MEMFREE(vtxdist);          Finley_setError(IO_ERROR, "saveVTK: spatial dimension 2 or 3 is supported only.");
113    MEMFREE(DOFNodes);          return;  
   MEMFREE(backwardBuffer);  
   MEMFREE(forwardBuffer);  
   
   if( myRank == 0)  
   {  
     char* tags = "</DataArray>\n</Points>\n<Cells>\n<DataArray Name=\"connectivity\" type=\"Int32\" " \  
                  "format=\"ascii\">\n" ;  
     MPI_File_iwrite_shared(fh,tags,strlen(tags),MPI_CHAR,&req);  
     MPI_Wait(&req,&status);  
114    }    }
115    MPIO_DEBUG(" Done Writing Coordinate Points ")    /*************************************************************************************/
116    
117    /* BEGIN CONNECTIVITY */    /* open the file and check handle */
118    
119    int NN = elements->ReferenceElement->Type->numNodes; /* num Nodes holding ref-element */    if (mpi_size > 1) {
120            #ifdef PASO_MPI
121              /* Collective Call */
122              #ifdef MPIO_HINTS
123                MPI_Info_create(&mpi_info);
124                /*  MPI_Info_set(mpi_info, "striping_unit",        "424288"); */
125                /*  MPI_Info_set(mpi_info, "striping_factor",      "16"); */
126                /*  MPI_Info_set(mpi_info, "collective_buffering", "true"); */
127                /*  MPI_Info_set(mpi_info, "cb_block_size",        "131072"); */
128                /*  MPI_Info_set(mpi_info, "cb_buffer_size",       "1048567"); */
129                /*  MPI_Info_set(mpi_info, "cb_nodes",             "8"); */
130                /*    MPI_Info_set(mpi_info, "access_style", "write_once, sequential"); */
131              
132                /*XFS only */
133                /*   MPI_Info_set(mpi_info, "direct_write",          "true"); */
134              #endif
135              ierr=MPI_File_open(mesh_p->Nodes->MPIInfo->comm, (char*)filename_p, amode,mpi_info, &mpi_fileHandle_p);
136              if (! ierr) {
137                  sprintf(error_msg, "saveVTK: File %s could not be opened for writing.", filename_p);
138                  Finley_setError(IO_ERROR,error_msg);
139              } else {
140                 MPI_File_set_view(mpi_fileHandle_p,MPI_DISPLACEMENT_CURRENT,MPI_CHAR, MPI_CHAR, "native" , mpi_info);
141              }
142            #endif
143      } else {
144            fileHandle_p = fopen(filename_p, "w");
145            if (fileHandle_p==NULL) {
146               sprintf(error_msg, "saveVTK: File %s could not be opened for writing.", filename_p);
147               Finley_setError(IO_ERROR,error_msg);
148             }
149      }
150      if (! Paso_MPIInfo_noError(mesh_p->Nodes->MPIInfo) ) return;
151      /*************************************************************************************/
152    
153    // Collective    /* find the mesh type to be written */
   MPIO_DEBUG(" Writing Connectivity... ")  
154    
155    size_t sz = numLocalCells*6*numVTKNodesPerElement + numLocalCells;    isCellCentered=TMPMEMALLOC(num_data,bool_t);
156    largebuf = MEMALLOC(sz,char);    max_len_names=0;
157    largebuf[0] = '\0';    if (!Finley_checkPtr(isCellCentered)) {
158    tsz=0;       for (i_data=0;i_data<num_data;++i_data) {
159    pos = 0;         if (! isEmpty(data_pp[i_data])) {
160    // numCells?           switch(getFunctionSpaceType(data_pp[i_data]) ) {
161    elementCache.values = MEMALLOC(numLocalCells,index_t);           case FINLEY_NODES:
162    if (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM)             nodetype = (nodetype == FINLEY_REDUCED_NODES) ? FINLEY_REDUCED_NODES : FINLEY_NODES;
163    {             isCellCentered[i_data]=FALSE;
164      for (i = 0; i < numCells; i++)             if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
165      {               elementtype=FINLEY_ELEMENTS;
166               } else {
167        if (elements->Id[i] >= elements->elementDistribution->vtxdist[i] &&  elements->Id[i] <= elements->elementDistribution->vtxdist[i+1] - 1 )               Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
168        {             }
169          for (j = 0; j < numVTKNodesPerElement; j++)             break;
170          {           case FINLEY_REDUCED_NODES:
171            sprintf(tmpbuf,"%d ",nodesGlobal[mesh_p->Nodes->toReduced[elements->Nodes[INDEX2(elements->ReferenceElement->Type->linearNodes[j], i, NN)]]]);             nodetype = FINLEY_REDUCED_NODES;
172            __STRCAT(largebuf,tmpbuf,tsz)             isCellCentered[i_data]=FALSE;
173               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
174                 elementtype=FINLEY_ELEMENTS;
175               } else {
176                 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
177               }
178               break;
179             case FINLEY_ELEMENTS:
180             case FINLEY_REDUCED_ELEMENTS:
181               nodetype = (nodetype == FINLEY_REDUCED_NODES) ? FINLEY_REDUCED_NODES :FINLEY_NODES;
182               isCellCentered[i_data]=TRUE;
183               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
184                 elementtype=FINLEY_ELEMENTS;
185               } else {
186                 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
187               }
188               break;
189             case FINLEY_FACE_ELEMENTS:
190             case FINLEY_REDUCED_FACE_ELEMENTS:
191               nodetype = (nodetype == FINLEY_REDUCED_NODES) ? FINLEY_REDUCED_NODES :FINLEY_NODES;
192               isCellCentered[i_data]=TRUE;
193               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_FACE_ELEMENTS) {
194                 elementtype=FINLEY_FACE_ELEMENTS;
195               } else {
196                 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
197               }
198               break;
199             case FINLEY_POINTS:
200               nodetype = (nodetype == FINLEY_REDUCED_NODES) ? FINLEY_REDUCED_NODES : FINLEY_NODES;
201               isCellCentered[i_data]=TRUE;
202               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_POINTS) {
203                 elementtype=FINLEY_POINTS;
204               } else {
205                 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
206               }
207               break;
208             case FINLEY_CONTACT_ELEMENTS_1:
209             case FINLEY_REDUCED_CONTACT_ELEMENTS_1:
210               nodetype = (nodetype == FINLEY_REDUCED_NODES) ? FINLEY_REDUCED_NODES :FINLEY_NODES;
211               isCellCentered[i_data]=TRUE;
212               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {
213                 elementtype=FINLEY_CONTACT_ELEMENTS_1;
214               } else {
215                 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
216               }
217               break;
218             case FINLEY_CONTACT_ELEMENTS_2:
219             case FINLEY_REDUCED_CONTACT_ELEMENTS_2:
220               nodetype = (nodetype == FINLEY_REDUCED_NODES) ? FINLEY_REDUCED_NODES :FINLEY_NODES;
221               isCellCentered[i_data]=TRUE;
222               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {
223                 elementtype=FINLEY_CONTACT_ELEMENTS_1;
224               } else {
225                 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
226               }
227               break;
228             default:
229               sprintf(error_msg,"saveVTK: unknown function space type %d",getFunctionSpaceType(data_pp[i_data]));
230               Finley_setError(TYPE_ERROR,error_msg);
231             }
232             if (isCellCentered[i_data]) {
233               write_celldata=TRUE;
234             } else {
235               write_pointdata=TRUE;
236             }
237             max_len_names =MAX(max_len_names,strlen(names_p[i_data]));
238           }
239         }
240      }
241      if (Finley_noError()) {
242    
243         /***************************************/
244    
245         /* select number of points and the mesh component */
246      
247         if (nodetype == FINLEY_REDUCED_NODES) {
248            myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
249            myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
250            globalNumPoints = Finley_NodeFile_getGlobalNumReducedNodes(mesh_p->Nodes);
251            globalNodeIndex= Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
252            nodeMapping=mesh_p->Nodes->reducedNodesMapping;
253         } else {
254            myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
255            myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
256            globalNumPoints = Finley_NodeFile_getGlobalNumNodes(mesh_p->Nodes);
257            globalNodeIndex= Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
258            nodeMapping=mesh_p->Nodes->nodesMapping;
259         }
260         myNumPoints = myLastNode - myFirstNode;
261         if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;
262         switch(elementtype) {
263           case FINLEY_ELEMENTS:
264              elements=mesh_p->Elements;
265              break;
266            case FINLEY_FACE_ELEMENTS:
267              elements=mesh_p->FaceElements;
268              break;
269            case FINLEY_POINTS:
270              elements=mesh_p->Points;
271              break;
272            case FINLEY_CONTACT_ELEMENTS_1:
273              elements=mesh_p->ContactElements;
274              break;
275         }
276         if (elements==NULL) {
277           Finley_setError(SYSTEM_ERROR,"saveVTK: undefined element file");
278         } else {
279           /* map finley element type to VTK element type */
280           numCells = elements->numElements;
281           globalNumCells = Finley_ElementFile_getGlobalNumElements(elements);
282           myNumCells= Finley_ElementFile_getMyNumElements(elements);
283           myFirstCell= Finley_ElementFile_getFirstElement(elements);
284           NN = elements->numNodes;
285           if (nodetype==FINLEY_REDUCED_NODES || nodetype==FINLEY_REDUCED_NODES) {
286              TypeId = elements->LinearReferenceElement->Type->TypeId;
287           } else {
288              TypeId = elements->ReferenceElement->Type->TypeId;
289           }
290           switch(TypeId) {
291            case Point1:
292            case Line2Face:
293            case Line3Face:
294            case Point1_Contact:
295            case Line2Face_Contact:
296            case Line3Face_Contact:
297              numCellFactor=1;
298              cellType = VTK_VERTEX;
299              numVTKNodesPerElement = 1;
300              strcpy(elemTypeStr, "VTK_VERTEX");
301              break;
302          
303            case Line2:
304            case Tri3Face:
305            case Rec4Face:
306            case Line2_Contact:
307            case Tri3_Contact:
308            case Tri3Face_Contact:
309            case Rec4Face_Contact:
310              numCellFactor=1;
311              cellType = VTK_LINE;
312              numVTKNodesPerElement = 2;
313              strcpy(elemTypeStr, "VTK_LINE");
314              break;
315          
316            case Tri3:
317            case Tet4Face:
318            case Tet4Face_Contact:
319              numCellFactor=1;
320              cellType = VTK_TRIANGLE;
321              numVTKNodesPerElement = 3;
322              strcpy(elemTypeStr, "VTK_TRIANGLE");
323              break;
324          
325            case Rec4:
326            case Hex8Face:
327            case Rec4_Contact:
328            case Hex8Face_Contact:
329              numCellFactor=1;
330              cellType = VTK_QUAD;
331              numVTKNodesPerElement = 4;
332              strcpy(elemTypeStr, "VTK_QUAD");
333              break;
334          
335            case Tet4:
336              numCellFactor=1;
337              cellType = VTK_TETRA;
338              numVTKNodesPerElement = 4;
339              strcpy(elemTypeStr, "VTK_TETRA");
340              break;
341          
342            case Hex8:
343              numCellFactor=1;
344              cellType = VTK_HEXAHEDRON;
345              numVTKNodesPerElement = 8;
346              strcpy(elemTypeStr, "VTK_HEXAHEDRON");
347              break;
348          
349            case Line3:
350            case Tri6Face:
351            case Rec8Face:
352            case Line3_Contact:
353            case Tri6Face_Contact:
354            case Rec8Face_Contact:
355              numCellFactor=1;
356              cellType = VTK_QUADRATIC_EDGE;
357              numVTKNodesPerElement = 3;
358              strcpy(elemTypeStr, "VTK_QUADRATIC_EDGE");
359              break;
360          
361            case Tri6:
362            case Tet10Face:
363            case Tri6_Contact:
364            case Tet10Face_Contact:
365              numCellFactor=1;
366              cellType = VTK_QUADRATIC_TRIANGLE;
367              numVTKNodesPerElement = 6;
368              strcpy(elemTypeStr, "VTK_QUADRATIC_TRIANGLE");
369              break;
370          
371            case Rec8:
372            case Hex20Face:
373            case Rec8_Contact:
374            case Hex20Face_Contact:
375              numCellFactor=1;
376              cellType = VTK_QUADRATIC_QUAD;
377              numVTKNodesPerElement = 8;
378              strcpy(elemTypeStr, "VTK_QUADRATIC_QUAD");
379              break;
380          
381            case Tet10:
382              numCellFactor=1;
383              cellType = VTK_QUADRATIC_TETRA;
384              numVTKNodesPerElement = 10;
385              strcpy(elemTypeStr, "VTK_QUADRATIC_TETRA");
386              break;
387          
388            case Hex20:
389              numCellFactor=1;
390              cellType = VTK_QUADRATIC_HEXAHEDRON;
391              numVTKNodesPerElement = 20;
392              strcpy(elemTypeStr, "VTK_QUADRATIC_HEXAHEDRON");
393              break;
394          
395            default:
396              sprintf(error_msg, "saveVTK: Element type %s is not supported by VTK",elements->ReferenceElement->Type->Name);
397              Finley_setError(VALUE_ERROR,error_msg);
398          }          }
399          __STRCAT(largebuf,newline,tsz)       }
         elementCache.values[pos++]=i;  
       }  
     }  
400    }    }
401    else if (VTK_QUADRATIC_HEXAHEDRON==cellType)    /***************************************/
   {  
     char tmpbuf2[20*20*2];  
   
     for (i = 0; i < numCells; i++)  
     {  
       if( elements->Id[i] >= elements->elementDistribution->vtxdist[myRank] && elements->Id[i] <= elements->elementDistribution->vtxdist[myRank+1]-1)  
       {  
         sprintf(tmpbuf2,"%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d\n",  
                 nodesGlobal[elements->Nodes[INDEX2(0, i, NN)]],  
                 nodesGlobal[elements->Nodes[INDEX2(1, i, NN)]],  
                 nodesGlobal[elements->Nodes[INDEX2(2, i, NN)]],  
                 nodesGlobal[elements->Nodes[INDEX2(3, i, NN)]],  
                 nodesGlobal[elements->Nodes[INDEX2(4, i, NN)]],  
                 nodesGlobal[elements->Nodes[INDEX2(5, i, NN)]],  
                 nodesGlobal[elements->Nodes[INDEX2(6, i, NN)]],  
                 nodesGlobal[elements->Nodes[INDEX2(7, i, NN)]],  
                 nodesGlobal[elements->Nodes[INDEX2(8, i, NN)]],  
                 nodesGlobal[elements->Nodes[INDEX2(9, i, NN)]],  
                 nodesGlobal[elements->Nodes[INDEX2(10, i, NN)]],  
                 nodesGlobal[elements->Nodes[INDEX2(11, i, NN)]],  
                 nodesGlobal[elements->Nodes[INDEX2(16, i, NN)]],  
                 nodesGlobal[elements->Nodes[INDEX2(17, i, NN)]],  
                 nodesGlobal[elements->Nodes[INDEX2(18, i, NN)]],  
                 nodesGlobal[elements->Nodes[INDEX2(19, i, NN)]],  
                 nodesGlobal[elements->Nodes[INDEX2(12, i, NN)]],  
                 nodesGlobal[elements->Nodes[INDEX2(13, i, NN)]],  
                 nodesGlobal[elements->Nodes[INDEX2(14, i, NN)]],  
                 nodesGlobal[elements->Nodes[INDEX2(15, i, NN)]]);  
         __STRCAT(largebuf,tmpbuf2,tsz)  
         elementCache.values[pos++]=i;  
       }  
     }  
   }  
   else if (numVTKNodesPerElement!=NN)  
   {  
402    
403      for (i = 0; i < numCells; i++)    /***************************************/
404      {    /*                                     */
405        if (elements->Id[i] >= elements->elementDistribution->vtxdist[i] &&  elements->Id[i] <= elements->elementDistribution->vtxdist[i+1] - 1 )    /*   allocate text buffer              */
406        {    /*                                     */
407          for (j = 0; j < numVTKNodesPerElement; j++)    len_txt_buffer= strlen(tags_header) + 3 * LEN_PRINTED_INT_FORMAT; /* header */
408          {    if (mpi_size > 1) len_txt_buffer=MAX(len_txt_buffer, myNumPoints * LEN_TMP_BUFFER);
409            sprintf(tmpbuf,"%d ", nodesGlobal[elements->Nodes[INDEX2(elements->ReferenceElement->Type->geoNodes[j], i, NN)]]);    if (mpi_size > 1) len_txt_buffer=MAX(len_txt_buffer, numCellFactor*myNumCells*(LEN_PRINTED_INT_FORMAT*numVTKNodesPerElement+1));
410            __STRCAT(largebuf,tmpbuf,tsz)    len_txt_buffer=MAX(len_txt_buffer,200+3*max_len_names);
411          }    len_txt_buffer=MAX(len_txt_buffer, strlen(tag_Float_DataArray) + LEN_PRINTED_INT_FORMAT + max_len_names);
412          __STRCAT(largebuf,newline,tsz)    if (mpi_size > 1) len_txt_buffer=MAX(len_txt_buffer, numCellFactor*myNumCells*LEN_PRINTED_FLOAT_TENSOR_FORMAT);
413          elementCache.values[pos++]=i;    if (mpi_size > 1) len_txt_buffer=MAX(len_txt_buffer, myNumPoints*LEN_PRINTED_FLOAT_TENSOR_FORMAT);
414        }    txt_buffer=TMPMEMALLOC(len_txt_buffer+1,char);
415      }    Finley_checkPtr(txt_buffer);
416    }    
417    else    if (Finley_noError()) {
   {  
     for(i = 0;i  < numCells ; i++)  
     {  
       if( elements->Id[i] >= elements->elementDistribution->vtxdist[myRank] && elements->Id[i] <= elements->elementDistribution->vtxdist[myRank+1]-1)  
       {  
         for (j = 0; j < numVTKNodesPerElement; j++)  
         {  
           sprintf(tmpbuf,"%d ", nodesGlobal[ elements->Nodes[INDEX2(j, i, NN) ] ] );  
           __STRCAT(largebuf,tmpbuf,tsz)  
         }  
         __STRCAT(largebuf,newline,tsz)  
         elementCache.values[pos++]=i;  
       }  
     }  
   }  
418    
419    elementCache.size = pos;       /* select number of points and the mesh component */
420    
421    largebuf[tsz] = '\0';       sprintf(txt_buffer,tags_header,globalNumPoints,numCellFactor*globalNumCells,3);
   MPI_File_write_ordered(fh,largebuf,tsz, MPI_CHAR, &status);  
   MEMFREE(largebuf);  
   MPIO_DEBUG(" Done Writing Connectivity ")  
   MPIO_DEBUG(" Writing Offsets & Types... ")  
   
   // Non-Collective  
   if( myRank == 0)  
   {  
     // write out the DataArray element for the offsets  
     char* tag1 = "</DataArray>\n<DataArray Name=\"offsets\" type=\"Int32\" format=\"ascii\">\n";  
     char* tag2 = "</DataArray>\n";  
     char *tag3 =  "<DataArray Name=\"types\" type=\"UInt8\" format=\"ascii\">\n";  
     char *tag4 = "</DataArray>\n</Cells>\n";  
   
     int n = numVTKNodesPerElement;  
   
     // allocate an upper bound on number of bytes needed    
     int sz=0;  
     int lg = log10(numGlobalCells * n) + 1;  
     sz += numGlobalCells*lg;  
     sz += numGlobalCells;    
     tsz = 0;  
   
     char* largebuf = MEMALLOC(sz  + strlen(tag1) + strlen(tag2) + strlen(tag3) + strlen(tag4),char);  
     largebuf[0] ='\0';  
     char tmp[10];  
     __STRCAT(largebuf,tag1,tsz)  
     for (i=numVTKNodesPerElement; i<=numGlobalCells*numVTKNodesPerElement; i+=numVTKNodesPerElement)  
     {  
       sprintf(tmp,"%d\n", i);  
       __STRCAT(largebuf,tmp,tsz)  
     }  
     __STRCAT(largebuf,tag2,tsz)  
     largebuf[tsz] = '\0';  
     MPI_File_iwrite_shared(fh,largebuf, tsz,MPI_CHAR,&req);  
     MPI_Wait(&req,&status);  
   
     // re-using buffer!  
     largebuf[0] = '\0';  
     tsz = 0;  
     __STRCAT(largebuf,tag3,tsz)  
     for (i=0; i<numGlobalCells; i++)  
     {  
       sprintf(tmp, "%d\n", cellType);  
       __STRCAT(largebuf,tmp,tsz)  
     }  
     __STRCAT(largebuf,tag4,tsz)  
     largebuf[tsz] = '\0';  
     MPI_File_iwrite_shared(fh,largebuf,tsz,MPI_CHAR,&req);  
     MPI_Wait(&req,&status);  
     MEMFREE(largebuf);  
   }  
422    
423    MPIO_DEBUG(" Done Writing Offsets & Types ")        if (mpi_size > 1) {
424              if ( my_mpi_rank == 0) {
425    // Write Cell data header Tags              #ifdef PASO_MPI
426    if(myRank == 0)                MPI_File_iwrite_shared(mpi_fileHandle_p,txt_buffer,strlen(txt_buffer),MPI_CHAR,&mpi_req);
427    {                MPI_Wait(&mpi_req,&mpi_status);
428      MPIO_DEBUG(" Writing Cell Data ...")              #endif
429      if( write_celldata)            }
430      {        } else {
431        char tmpBuf[80];           fprintf(fileHandle_p,txt_buffer);
432        char header[600];        }
433        // mark the active data arrays  
434        bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;        /* write the nodes */
435        sprintf(tmpBuf, "<CellData");        
436        strcat(header,tmpBuf);        if (mpi_size > 1) {
437        for (i_data =0 ;i_data<num_data;++i_data)  
438        {           txt_buffer[0] = '\0';
439          if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data])           txt_buffer_in_use=0;
440          {           if (nDim==2) {
441            // if the rank == 0:   --> scalar data              for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
442            // if the rank == 1:   --> vector data                 if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
443            // if the rank == 2:   --> tensor data                   sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,
444                                        mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
445                                        mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
446                                        0.);
447                     __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use);
448                   }
449                }      
450             } else {
451                for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
452                   if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
453                     sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,
454                                                      mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
455                                                      mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
456                                                      mesh_p->Nodes->Coordinates[INDEX2(2, i, nDim)]);
457                     __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use);
458                   }
459                }    
460      
461             }
462             #ifdef PASO_MPI
463                MPI_File_write_ordered(mpi_fileHandle_p, txt_buffer,txt_buffer_in_use, MPI_CHAR, &mpi_status);
464             #endif    
465          } else {
466             if (nDim==2) {
467                for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
468                   if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
469                     fprintf(fileHandle_p,FLOAT_VECTOR_FORMAT,
470                                          mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
471                                          mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
472                                          0.);
473                   }
474                }      
475             } else {
476                for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
477                   if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
478                     fprintf(fileHandle_p,FLOAT_VECTOR_FORMAT,
479                                                  mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
480                                                  mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
481                                                  mesh_p->Nodes->Coordinates[INDEX2(2, i, nDim)]);
482                   }
483                }    
484      
485             }
486          }
487    
488          /* close the Points and open connectivity */
489    
490          if (mpi_size > 1) {
491              if ( my_mpi_rank == 0) {
492                 #ifdef PASO_MPI
493                    MPI_File_iwrite_shared(mpi_fileHandle_p, tags_End_Points_and_Start_Conn, strlen(tags_End_Points_and_Start_Conn), MPI_CHAR, &mpi_req);
494                    MPI_Wait(&mpi_req,&mpi_status);
495                 #endif
496              } else {
497                 fprintf(fileHandle_p,tags_End_Points_and_Start_Conn);
498              }
499        
500           }
501    
502         /* write the cells */
503         if (nodetype == FINLEY_REDUCED_NODES) {
504            node_index=elements->ReferenceElement->Type->linearNodes;
505         } else if (VTK_QUADRATIC_HEXAHEDRON==cellType) {
506            node_index=VTK_QUADRATIC_HEXAHEDRON_INDEX;
507         } else if (numVTKNodesPerElement!=NN) {
508            node_index=elements->ReferenceElement->Type->geoNodes;
509         } else {
510            node_index=NULL;
511         }
512    
513         if ( mpi_size > 1) {
514            txt_buffer[0] = '\0';
515            txt_buffer_in_use=0;
516            if (node_index == NULL) {
517               for (i = 0; i < numCells; i++) {
518                  if (elements->Owner[i] == my_mpi_rank) {
519                     for (j = 0; j < numVTKNodesPerElement; j++) {
520                         sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(j, i, NN)]]);
521                         __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use)
522                     }
523                     __STRCAT(txt_buffer,NEWLINE,txt_buffer_in_use)
524                  }
525               }
526            } else {
527               for (i = 0; i < numCells; i++) {
528                  if (elements->Owner[i] == my_mpi_rank) {
529                     for (j = 0; j < numVTKNodesPerElement; j++) {
530                         sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(node_index[j], i, NN)]]);
531                         __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use)
532                     }
533                     __STRCAT(txt_buffer,NEWLINE,txt_buffer_in_use)
534                  }
535               }
536            }
537            #ifdef PASO_MPI
538               MPI_File_write_ordered(mpi_fileHandle_p,txt_buffer,txt_buffer_in_use, MPI_CHAR, &mpi_status);
539            #endif    
540         } else {
541            if (node_index == NULL) {
542               for (i = 0; i < numCells; i++) {
543                  for (j = 0; j < numVTKNodesPerElement; j++) {
544                     fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(j, i, NN)]]);
545                   }
546                  fprintf(fileHandle_p,NEWLINE);
547               }
548            } else {
549               for (i = 0; i < numCells; i++) {
550                  for (j = 0; j < numVTKNodesPerElement; j++) {
551                     fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(node_index[j], i, NN)]]);
552                   }
553                  fprintf(fileHandle_p,NEWLINE);
554               }
555            }
556    
557         }
558        
559         /* finalize the connection and start the offset section */
560         if (mpi_size > 1) {
561            if( my_mpi_rank == 0) {
562               #ifdef PASO_MPI
563                  MPI_File_iwrite_shared(mpi_fileHandle_p,tags_End_Conn_and_Start_Offset,strlen(tags_End_Conn_and_Start_Offset),MPI_CHAR,&mpi_req);
564                  MPI_Wait(&mpi_req,&mpi_status);
565               #endif
566            }
567         } else {
568            fprintf(fileHandle_p,tags_End_Conn_and_Start_Offset);
569         }
570    
571        /* write the offsets */
572          
573         if ( mpi_size > 1) {
574            txt_buffer[0] = '\0';
575            txt_buffer_in_use=0;
576            for (i=numVTKNodesPerElement*myFirstCell*numCellFactor; i<=numCells*numVTKNodesPerElement*numCellFactor; i+=numVTKNodesPerElement) {
577               sprintf(tmp_buffer, INT_NEWLINE_FORMAT, i);
578               __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use);
579             }
580             #ifdef PASO_MPI
581                MPI_File_write_ordered(mpi_fileHandle_p,txt_buffer,txt_buffer_in_use, MPI_CHAR, &mpi_status);
582             #endif    
583         } else {
584            for (i=0; i<=numCells*numVTKNodesPerElement*numCellFactor; i+=numVTKNodesPerElement) {
585               fprintf(fileHandle_p, INT_NEWLINE_FORMAT, i);
586            }
587        
588         }
589         /* finalize the offset section and start the type section */
590         if ( mpi_size > 1) {
591            if ( my_mpi_rank == 0) {
592               #ifdef PASO_MPI
593                  MPI_File_iwrite_shared(mpi_fileHandle_p,tags_End_Offset_and_Start_Type,strlen(tags_End_Offset_and_Start_Type),MPI_CHAR,&mpi_req);
594                  MPI_Wait(&mpi_req,&mpi_status);
595               #endif
596            }
597        } else {
598           fprintf(fileHandle_p,tags_End_Offset_and_Start_Type);
599        }
600        
601          
602         /* write element type */
603         sprintf(tmp_buffer, INT_NEWLINE_FORMAT, cellType);
604         if ( mpi_size > 1) {
605            txt_buffer[0] = '\0';
606            txt_buffer_in_use=0;
607            for (i=0; i<=numCells*numCellFactor; i++) __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use);
608             #ifdef PASO_MPI
609                MPI_File_write_ordered(mpi_fileHandle_p,txt_buffer,txt_buffer_in_use, MPI_CHAR, &mpi_status);
610             #endif    
611         } else {
612            for (i=0; i<=numCells*numCellFactor; i++) fprintf(fileHandle_p, tmp_buffer);
613         }
614         /* finalize cell information */
615         if ( mpi_size > 1) {
616            if ( my_mpi_rank == 0) {
617               #ifdef PASO_MPI
618                  MPI_File_iwrite_shared(mpi_fileHandle_p,tags_End_Type_And_Cells,strlen(tags_End_Type_And_Cells),MPI_CHAR,&mpi_req);
619                  MPI_Wait(&mpi_req,&mpi_status);
620               #endif
621            }
622        } else {
623           fprintf(fileHandle_p,tags_End_Type_And_Cells);
624        }
625     }
626    
627     /* Write cell data */
628     if (write_celldata && Finley_noError()) {
629          /* mark the active data arrays */
630          set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
631          strcat(txt_buffer, "<CellData");
632          for (i_data =0 ;i_data<num_data;++i_data) {
633            if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data]) {
634              /* if the rank == 0:   --> scalar data */
635              /* if the rank == 1:   --> vector data */
636              /* if the rank == 2:   --> tensor data */
637    
638            switch(getDataPointRank(data_pp[i_data]))            switch(getDataPointRank(data_pp[i_data])) {
           {  
639            case 0:            case 0:
640              if (! set_scalar)              if (! set_scalar) {
641              {                strcat(txt_buffer," Scalars=\"");
642                sprintf(tmpBuf," Scalars=\"%s\"",names_p[i_data]);                strcat(txt_buffer,names_p[i_data]);
643                strcat(header,tmpBuf);                strcat(txt_buffer,"\"");
644                set_scalar=TRUE;                set_scalar=TRUE;
645              }              }
646              break;              break;
647            case 1:            case 1:
648              if (! set_vector)              if (! set_vector) {
649              {                strcat(txt_buffer," Vectors=\"");
650                sprintf(tmpBuf," Vectors=\"%s\"",names_p[i_data]);                strcat(txt_buffer,names_p[i_data]);
651            strcat(header,tmpBuf);                strcat(txt_buffer,"\"");
652                set_vector=TRUE;                set_vector=TRUE;
653              }              }
654              break;              break;
655            case 2:            case 2:
656              if (! set_tensor)              if (! set_tensor) {
657              {                strcat(txt_buffer," Tensors=\"");
658                sprintf(tmpBuf," Tensors=\"%s\"",names_p[i_data]);                strcat(txt_buffer,names_p[i_data]);
659            strcat(header,tmpBuf);                strcat(txt_buffer,"\"");
660                set_tensor=TRUE;                set_tensor=TRUE;
661              }              }
662              break;              break;
# Line 804  void Finley_Mesh_saveVTK_MPIO(const char Line 667  void Finley_Mesh_saveVTK_MPIO(const char
667            }            }
668          }          }
669        }        }
670        strcat(header, ">\n");        strcat(txt_buffer, ">\n");
671        MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);        if ( mpi_size > 1) {
672        MPI_Wait(&req,&status);          if ( my_mpi_rank == 0) {
673      }             #ifdef PASO_MPI
674    }                MPI_File_iwrite_shared(mpi_fileHandle_p,txt_buffer,strlen(txt_buffer),MPI_CHAR,&mpi_req);
675                  MPI_Wait(&mpi_req,&mpi_status);
676    // write actual data (collective)             #endif
677    if(write_celldata)          }
678    {        } else {
679      for (i_data =0 ;i_data<num_data;++i_data)            fprintf(fileHandle_p,txt_buffer);
680      {        }
681        if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data])        /* write the arrays */
682        {        for (i_data =0 ;i_data<num_data;++i_data) {
683          numPointsPerSample = elements->ReferenceElement->numQuadNodes;           if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data]) {
684          rank = getDataPointRank(data_pp[i_data]);              txt_buffer[0] = '\0';
685          nComp = getDataPointSize(data_pp[i_data]);              txt_buffer_in_use=0;
686          nCompReqd=1;   // the number of components required by vtk              numPointsPerSample=getNumDataPointsPerSample(data_pp[i_data]);
687          shape=0;              rank = getDataPointRank(data_pp[i_data]);
688          if (rank == 0)              nComp = getDataPointSize(data_pp[i_data]);
689          {              nCompReqd=1;   /* the number of components mpi_required by vtk */
690            nCompReqd = 1;              shape=0;
691          }              if (rank == 0) {
692          else if (rank == 1)                nCompReqd = 1;
693          {              } else if (rank == 1) {
694            shape=getDataPointShape(data_pp[i_data], 0);                shape=getDataPointShape(data_pp[i_data], 0);
695            if  (shape>3)                if  (shape>3) {
696            {                  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");  
             return;  
           }  
           nCompReqd = 3;  
         }  
         else  
         {  
           shape=getDataPointShape(data_pp[i_data], 0);  
           if  (shape>3 || shape != getDataPointShape(data_pp[i_data], 1))  
           {  
             Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");  
             return;  
           }  
           nCompReqd = 9;  
         }  
   
         if( myRank == 0)  
         {  
           char header[250];  
           sprintf(header,"<DataArray Name=\"%s\" type=\"Float64\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);  
           MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);  
           MPI_Wait(&req,&status);  
         }  
   
         // Write the actual data  
         char tmpbuf[15];  
         char* largebuf = MEMALLOC(nCompReqd*numLocalCells*15 + numLocalCells*nCompReqd + nCompReqd + 15,char);  
         largebuf[0] = '\0';  
         size_t tsz = 0;  
   
         double sampleAvg[nComp];  
   
         for (k=0; k<elementCache.size; k++)  
         {  
           i = elementCache.values[k];  
   
           values = getSampleData(data_pp[i_data], i);  
           // averaging over the number of points in the sample  
           for (n=0; n<nComp; n++)  
           {  
             if (isExpanded(data_pp[i_data])) {  
                rtmp = 0.;  
                for (j=0; j<numPointsPerSample; j++) rtmp += values[INDEX2(n,j,nComp)];  
                sampleAvg[n] = rtmp/numPointsPerSample;  
             } else {  
                sampleAvg[n] = values[n];  
             }  
           }  
           // if the number of required components is more than the number  
           // of actual components, pad with zeros  
   
           // probably only need to get shape of first element  
           // write the data different ways for scalar, vector and tensor  
           if (nCompReqd == 1)  
           {  
             sprintf(tmpbuf, " %e", sampleAvg[0]);  
             __STRCAT(largebuf,tmpbuf,tsz)  
           }  
           else if (nCompReqd == 3)  
           {  
             // write out the data  
             for (m=0; m<shape; m++)  
             {  
               sprintf(tmpbuf, " %e", sampleAvg[m]);  
               __STRCAT(largebuf,tmpbuf,tsz)  
             }  
             for (m=0; m<nCompReqd-shape; m++)  
             {  
               __STRCAT(largebuf,zero,tsz)  
             }  
           }  
           else if (nCompReqd == 9)  
           {  
             // tensor data, so have a 3x3 matrix to output as a row  
             // of 9 data points  
             count = 0;  
             for (m=0; m<shape; m++)  
             {  
               for (n=0; n<shape; n++)  
               {  
                 sprintf(tmpbuf, " %e", sampleAvg[count]);  
                 __STRCAT(largebuf,tmpbuf,tsz)  
                 count++;  
               }  
               for (n=0; n<3-shape; n++)  
               {  
                 __STRCAT(largebuf,zero,tsz)  
               }  
             }  
             for (m=0; m<3-shape; m++)  
               for (n=0; n<3; n++)  
               {  
                 __STRCAT(largebuf,zero,tsz)  
697                }                }
698            }                nCompReqd = 3;
699            __STRCAT(largebuf,newline,tsz)              } else {
700          }                shape=getDataPointShape(data_pp[i_data], 0);
701          largebuf[tsz] = '\0';                if  (shape>3 || shape != getDataPointShape(data_pp[i_data], 1)) {
702          MPI_File_write_ordered(fh,largebuf,tsz,MPI_CHAR,&status);                  Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
703          MEMFREE(largebuf);                }
704          if( myRank == 0)                nCompReqd = 9;
705          {              }
706            char *tag = "</DataArray>\n";              if (Finley_noError()) {
707            MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);                 sprintf(txt_buffer,tag_Float_DataArray,names_p[i_data], nCompReqd);
708            MPI_Wait(&req,&status);                 if ( mpi_size > 1) {
709                     if ( my_mpi_rank == 0) {
710                        #ifdef PASO_MPI
711                           MPI_File_iwrite_shared(mpi_fileHandle_p,txt_buffer,strlen(txt_buffer),MPI_CHAR,&mpi_req);
712                           MPI_Wait(&mpi_req,&mpi_status);
713                        #endif
714                     }
715                   } else {
716                       fprintf(fileHandle_p,txt_buffer);
717                   }
718                   for (i=0; i<numCells; i++) {
719                       if (elements->Owner[i] == my_mpi_rank) {
720                          values = getSampleData(data_pp[i_data], i);
721                          /* averaging over the number of points in the sample */
722                          for (k=0; k<MIN(nComp,NCOMP_MAX); k++) {
723                             if (isExpanded(data_pp[i_data])) {
724                               rtmp = 0.;
725                               for (j=0; j<numPointsPerSample; j++) rtmp += values[INDEX2(k,j,nComp)];
726                               sampleAvg[k] = rtmp/numPointsPerSample;
727                            } else {
728                               sampleAvg[k] = values[k];
729                            }
730                          }
731                          /* if the number of mpi_required components is more than the number
732                          * of actual components, pad with zeros
733                          */
734                          /* probably only need to get shape of first element */
735                          /* write the data different ways for scalar, vector and tensor */
736                          if (nCompReqd == 1) {
737                            sprintf(tmp_buffer,FLOAT_SCALAR_FORMAT,sampleAvg[0]);
738                          } else if (nCompReqd == 3) {
739                            if (shape==1) {
740                             sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,sampleAvg[0],0.,0.);
741                            } else if (shape==2) {
742                             sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,sampleAvg[0],sampleAvg[1],0.);
743                            } else if (shape==3) {
744                             sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,sampleAvg[0],sampleAvg[2],sampleAvg[3]);
745                            }
746                          } else if (nCompReqd == 9) {
747                            if (shape==1) {
748                             sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,sampleAvg[0],0.,0.,
749                                                                    0.,0.,0.,
750                                                                    0.,0.,0.);
751                            } else if (shape==2) {
752                             sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,sampleAvg[0],sampleAvg[1],0.,
753                                                                    sampleAvg[2],sampleAvg[3],0.,
754                                                                    0.,0.,0.);
755                            } else if (shape==3) {
756                             sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,sampleAvg[0],sampleAvg[1],sampleAvg[2],
757                                                                    sampleAvg[3],sampleAvg[4],sampleAvg[5],
758                                                                    sampleAvg[6],sampleAvg[7],sampleAvg[8]);
759                            }
760                          }
761                          if ( mpi_size > 1) {
762                            __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use);
763                          } else {
764                            fprintf(fileHandle_p,tmp_buffer);
765                          }
766                      }
767                   }
768                   if ( mpi_size > 1) {
769                         #ifdef PASO_MPI
770                            MPI_File_write_ordered(mpi_fileHandle_p,txt_buffer,txt_buffer_in_use, MPI_CHAR, &mpi_status);
771                         #endif    
772                         if ( my_mpi_rank == 0) {
773                            #ifdef PASO_MPI
774                               MPI_File_iwrite_shared(mpi_fileHandle_p,tag_End_DataArray,strlen(tag_End_DataArray),MPI_CHAR,&mpi_req);
775                               MPI_Wait(&mpi_req,&mpi_status);
776                            #endif
777                         }
778                   } else {
779                       fprintf(fileHandle_p,tag_End_DataArray);
780                   }
781                }
782             }
783          }
784          if ( mpi_size > 1) {
785            if ( my_mpi_rank == 0) {
786               #ifdef PASO_MPI
787                  MPI_File_iwrite_shared(mpi_fileHandle_p,tag_End_CellData,strlen(tag_End_CellData),MPI_CHAR,&mpi_req);
788                  MPI_Wait(&mpi_req,&mpi_status);
789               #endif
790          }          }
791          } else {
792              fprintf(fileHandle_p,tag_End_CellData);
793        }        }
     }  
     // closing celldata tag  
     if(myRank == 0)  
     {  
       char* tag =  "</CellData>\n";  
       MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);  
       MPI_Wait(&req,&status);  
     }  
   
     MPIO_DEBUG(" Done Writing Cell Data ")  
794    }    }
795      /* point data */
796      if (write_pointdata && Finley_noError()) {
797          /* mark the active data arrays */
798          set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
799          strcat(txt_buffer, "<PointData");
800          for (i_data =0 ;i_data<num_data;++i_data) {
801            if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data]) {
802              /* if the rank == 0:   --> scalar data */
803              /* if the rank == 1:   --> vector data */
804              /* if the rank == 2:   --> tensor data */
805    
806              switch(getDataPointRank(data_pp[i_data])) {
   // Write Point Data Header Tags  
   if( myRank == 0)  
   {  
     char header[600];  
     char tmpBuf[50];  
   
     if (write_pointdata)  
     {  
       MPIO_DEBUG(" Writing Pointdata... ")  
       // mark the active data arrays  
       bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;  
       sprintf(header, "<PointData");  
       for (i_data =0 ;i_data<num_data;++i_data)  
       {  
         if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])  
         {  
           // if the rank == 0:   --> scalar data  
           // if the rank == 1:   --> vector data  
           // if the rank == 2:   --> tensor data  
   
           switch(getDataPointRank(data_pp[i_data]))  
           {  
807            case 0:            case 0:
808              if (! set_scalar)              if (! set_scalar) {
809              {                strcat(txt_buffer," Scalars=\"");
810                sprintf(tmpBuf," Scalars=\"%s\"",names_p[i_data]);                strcat(txt_buffer,names_p[i_data]);
811                strcat(header,tmpBuf);                strcat(txt_buffer,"\"");
812                set_scalar=TRUE;                set_scalar=TRUE;
813              }              }
814              break;              break;
815            case 1:            case 1:
816              if (! set_vector)              if (! set_vector) {
817              {                strcat(txt_buffer," Vectors=\"");
818                sprintf(tmpBuf," Vectors=\"%s\"",names_p[i_data]);                strcat(txt_buffer,names_p[i_data]);
819                strcat(header,tmpBuf);                strcat(txt_buffer,"\"");
820                set_vector=TRUE;                set_vector=TRUE;
821              }              }
822              break;              break;
823            case 2:            case 2:
824              if (! set_tensor)              if (! set_tensor) {
825              {                strcat(txt_buffer," Tensors=\"");
826                sprintf(tmpBuf," Tensors=\"%s\"",names_p[i_data]);                strcat(txt_buffer,names_p[i_data]);
827                strcat(header,tmpBuf);                strcat(txt_buffer,"\"");
828                set_tensor=TRUE;                set_tensor=TRUE;
829              }              }
830              break;              break;
# Line 1005  void Finley_Mesh_saveVTK_MPIO(const char Line 835  void Finley_Mesh_saveVTK_MPIO(const char
835            }            }
836          }          }
837        }        }
838        strcat(header, ">\n");        strcat(txt_buffer, ">\n");
839        MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);        if ( mpi_size > 1) {
840        MPI_Wait(&req,&status);          if ( my_mpi_rank == 0) {
841      }             #ifdef PASO_MPI
842    }                MPI_File_iwrite_shared(mpi_fileHandle_p,txt_buffer,strlen(txt_buffer),MPI_CHAR,&mpi_req);
843                  MPI_Wait(&mpi_req,&mpi_status);
844    // write actual data             #endif
845    if(write_pointdata)          }
846    {        } else {
847      for (i_data =0 ;i_data<num_data;++i_data)            fprintf(fileHandle_p,txt_buffer);
848      {        }
849        if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])        /* write the arrays */
850        {        for (i_data =0 ;i_data<num_data;++i_data) {
851          if (Finley_Assemble_reducedIntegrationOrder(data_pp[i_data])) {           if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data]) {
852             numPointsPerSample=elements->ReferenceElementReducedOrder->numQuadNodes;              txt_buffer[0] = '\0';
853          } else {              txt_buffer_in_use=0;
854             numPointsPerSample=elements->ReferenceElement->numQuadNodes;              numPointsPerSample=getNumDataPointsPerSample(data_pp[i_data]);
855          }              rank = getDataPointRank(data_pp[i_data]);
856          rank = getDataPointRank(data_pp[i_data]);              nComp = getDataPointSize(data_pp[i_data]);
857          nComp = getDataPointSize(data_pp[i_data]);              nCompReqd=1;   /* the number of components mpi_required by vtk */
858          nCompReqd=1;   // the number of components required by vtk              shape=0;
859          shape=0;              if (rank == 0) {
860          if (rank == 0)                nCompReqd = 1;
861          {              } else if (rank == 1) {
862            nCompReqd = 1;                shape=getDataPointShape(data_pp[i_data], 0);
863          }                if  (shape>3) {
864          else if (rank == 1)                  Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");
         {  
           shape=getDataPointShape(data_pp[i_data], 0);  
           if  (shape>3)  
           {  
             Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");  
             return;  
           }  
           nCompReqd = 3;  
         }  
         else  
         {  
           shape=getDataPointShape(data_pp[i_data], 0);  
           if  (shape>3 || shape != getDataPointShape(data_pp[i_data], 1))  
           {  
             Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");  
             return;  
           }  
           nCompReqd = 9;  
         }  
   
         if( myRank == 0)  
         {  
           char header[250];  
           sprintf(header, "<DataArray Name=\"%s\" type=\"Float64\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);  
           MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);  
           MPI_Wait(&req,&status);  
         }  
         // write out the data  
         // if the number of required components is more than the number  
         // of actual components, pad with zeros  
   
         char tmpbuf[15];  
         char* largebuf = MEMALLOC(nCompReqd*numLocalNodes*15 + numLocal*nCompReqd + nCompReqd + 15,char);  
         largebuf[0] = '\0';  
         bool_t do_write=TRUE;  
         size_t tsz = 0;  
   
         for(k=0;k < nodeCache.size;k++)  
         {  
           i = nodeCache.values[k];  
   
           if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)  
           {  
             if (mesh_p->Nodes->toReduced[i]>=0)  
             {  
               switch(getFunctionSpaceType(data_pp[i_data]))  
               {  
               case FINLEY_DEGREES_OF_FREEDOM:  
                 values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);  
                 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;  
             }  
           }  
           else  
           {  
             do_write=TRUE;  
             switch(getFunctionSpaceType(data_pp[i_data]))  
             {  
             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;  
             }  
           }  
           // write the data different ways for scalar, vector and tensor  
           if (do_write)  
           {  
             if (nCompReqd == 1)  
             {  
               sprintf(tmpbuf," %e", values[0]);  
               __STRCAT(largebuf,tmpbuf,tsz)  
             }  
             else if (nCompReqd == 3)  
             {  
               for (m=0; m<shape; m++)  
               {  
   
                 sprintf(tmpbuf," %e",values[m]);  
                 __STRCAT(largebuf,tmpbuf,tsz)  
865                }                }
866                for (m=0; m<nCompReqd-shape; m++)                nCompReqd = 3;
               {  
                 __STRCAT(largebuf,zero,tsz)  
               }  
             }  
             else if (nCompReqd == 9)  
             {  
               // tensor data, so have a 3x3 matrix to output as a row  
               //  of 9 data points  
               count = 0;  
               for (m=0; m<shape; m++)  
               {  
                 for (n=0; n<shape; n++)  
                 {  
                   sprintf(tmpbuf," %e", values[count]);  
                   __STRCAT(largebuf,tmpbuf,tsz)  
                   count++;  
                 }  
                 for (n=0; n<3-shape; n++)  
                 {  
                   __STRCAT(largebuf,zero,tsz)  
                 }  
               }  
               for (m=0; m<3-shape; m++)  
               {  
                 for (n=0; n<3; n++)  
                 {  
                   __STRCAT(largebuf,zero,tsz)  
                 }  
               }  
             }  
             __STRCAT(largebuf,newline,tsz)  
           }  
   
         }  
         // Write out local data  
   
         largebuf[tsz] = '\0';  
         MPI_File_write_ordered(fh,largebuf,tsz,MPI_CHAR,&status);  
         MEMFREE(largebuf);  
         if( myRank == 0)  
         {  
           char *tag = "</DataArray>\n";  
           MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);  
           MPI_Wait(&req,&status);  
         }  
       }  
     }  
     // Finish off with closing tag  
     if(myRank == 0)  
     {  
       char* tag =  "</PointData>\n";  
       MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);  
       MPI_Wait(&req,&status);  
     }  
     MPIO_DEBUG(" Done Writing Pointdata ")  
   }  
   // end write_pointdata  
   
   // tag and bag...    
   if (myRank == 0)  
   {  
     char *footer = "</Piece>\n</UnstructuredGrid>\n</VTKFile>";  
     MPI_File_iwrite_shared(fh,footer,strlen(footer),MPI_CHAR,&req);  
     MPI_Wait(&req,&status);  
   }  
   
   MEMFREE(nodesGlobal);  
   MEMFREE(nodeCache.values);  
   MEMFREE(elementCache.values);  
 #ifdef MPIO_HINTS  
   MPI_Info_free(&infoHints);  
 #undef MPIO_HINTS  
 #endif  
   MPI_File_close(&fh);  
   MPIO_DEBUG(" ***** Exit saveVTK ***** ")  
 #undef __STRCAT  
 }  
   
 #undef MPIO_DEBUG  
 #else  
   
   
   
   
 void Finley_Mesh_saveVTK(const char * filename_p, Finley_Mesh *mesh_p, const dim_t num_data,char* *names_p, escriptDataC* *data_pp)  
 {  
   #define NCOMP_MAX 9  
   char error_msg[LenErrorMsg_MAX];  
   double sampleAvg[NCOMP_MAX];  
   /* if there is no mesh we just return */  
   
   int i, j, k, numVTKNodesPerElement,i_data,m, count, n, rank,shape, numPoints, cellType, numCells,  
   nDim, numPointsPerSample, nComp, nCompReqd, NN;  
   index_t j2;  
   int nodetype=FINLEY_DEGREES_OF_FREEDOM;  
   int elementtype=FINLEY_UNKNOWN;  
   double* values, rtmp;  
   char elemTypeStr[32];  
   FILE * fileHandle_p = NULL;  
   bool_t do_write, *isCellCentered=NULL,write_celldata=FALSE,write_pointdata=FALSE;  
   Finley_ElementFile* elements=NULL;  
   ElementTypeId TypeId;  
   
   /* open the file and check handle */  
   if (mesh_p==NULL) return;  
   
   fileHandle_p = fopen(filename_p, "w");  
   if (fileHandle_p==NULL)  
   {  
     sprintf(error_msg, "saveVTK: File %s could not be opened for writing.", filename_p);  
     Finley_setError(IO_ERROR,error_msg);  
     return;  
   }  
   /* find the mesh type to be written */  
   isCellCentered=TMPMEMALLOC(num_data,bool_t);  
   
   
   if (Finley_checkPtr(isCellCentered)) {  
      fclose(fileHandle_p);  
      return;  
   }  
   for (i_data=0;i_data<num_data;++i_data)  
   {  
     if (! isEmpty(data_pp[i_data]))  
     {  
       switch(getFunctionSpaceType(data_pp[i_data]))  
       {  
       case FINLEY_DEGREES_OF_FREEDOM:  
         nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;  
         if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)  
         {  
           elementtype=FINLEY_ELEMENTS;  
         }  
         else  
         {  
           Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");  
           fclose(fileHandle_p);  
           TMPMEMFREE(isCellCentered);  
           return;  
         }  
         isCellCentered[i_data]=FALSE;  
         break;  
       case FINLEY_REDUCED_DEGREES_OF_FREEDOM:  
         nodetype = FINLEY_REDUCED_DEGREES_OF_FREEDOM;  
         if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)  
         {  
           elementtype=FINLEY_ELEMENTS;  
         }  
         else  
         {  
           Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");  
           fclose(fileHandle_p);  
           TMPMEMFREE(isCellCentered);  
           return;  
         }  
         isCellCentered[i_data]=FALSE;  
         break;  
       case FINLEY_NODES:  
         nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;  
         if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)  
         {  
           elementtype=FINLEY_ELEMENTS;  
         }  
         else  
         {  
           Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");  
           fclose(fileHandle_p);  
           TMPMEMFREE(isCellCentered);  
           return;  
         }  
         isCellCentered[i_data]=FALSE;  
         break;  
       case FINLEY_ELEMENTS:  
       case FINLEY_REDUCED_ELEMENTS:  
         nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;  
         if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)  
         {  
           elementtype=FINLEY_ELEMENTS;  
         }  
         else  
         {  
           Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");  
           fclose(fileHandle_p);  
           return;  
           TMPMEMFREE(isCellCentered);  
         }  
         isCellCentered[i_data]=TRUE;  
         break;  
       case FINLEY_FACE_ELEMENTS:  
       case FINLEY_REDUCED_FACE_ELEMENTS:  
         nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;  
         if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_FACE_ELEMENTS)  
         {  
           elementtype=FINLEY_FACE_ELEMENTS;  
         }  
         else  
         {  
           Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");  
           fclose(fileHandle_p);  
           TMPMEMFREE(isCellCentered);  
           return;  
         }  
         isCellCentered[i_data]=TRUE;  
         break;  
       case FINLEY_POINTS:  
         nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;  
         if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_POINTS)  
         {  
           elementtype=FINLEY_POINTS;  
         }  
         else  
         {  
           Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");  
           fclose(fileHandle_p);  
           TMPMEMFREE(isCellCentered);  
           return;  
         }  
         isCellCentered[i_data]=TRUE;  
         break;  
       case FINLEY_CONTACT_ELEMENTS_1:  
       case FINLEY_REDUCED_CONTACT_ELEMENTS_1:  
         nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;  
         if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1)  
         {  
           elementtype=FINLEY_CONTACT_ELEMENTS_1;  
         }  
         else  
         {  
           Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");  
           fclose(fileHandle_p);  
           TMPMEMFREE(isCellCentered);  
           return;  
         }  
         isCellCentered[i_data]=TRUE;  
         break;  
       case FINLEY_CONTACT_ELEMENTS_2:  
       case FINLEY_REDUCED_CONTACT_ELEMENTS_2:  
         nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;  
         if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1)  
         {  
           elementtype=FINLEY_CONTACT_ELEMENTS_1;  
         }  
         else  
         {  
           Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");  
           fclose(fileHandle_p);  
           TMPMEMFREE(isCellCentered);  
           return;  
         }  
         isCellCentered[i_data]=TRUE;  
         break;  
       default:  
         sprintf(error_msg,"saveVTK: Finley does not know anything about function space type %d",getFunctionSpaceType(data_pp[i_data]));  
         Finley_setError(TYPE_ERROR,error_msg);  
         fclose(fileHandle_p);  
         TMPMEMFREE(isCellCentered);  
         return;  
       }  
       if (isCellCentered[i_data])  
       {  
         write_celldata=TRUE;  
       }  
       else  
       {  
         write_pointdata=TRUE;  
       }  
     }  
   }  
   /* select nomber of points and the mesh component */  
   numPoints = mesh_p->Nodes->numNodes;  
   if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)  
   {  
     numPoints = mesh_p->Nodes->reducedNumNodes;  
   }  
   else  
   {  
     numPoints = mesh_p->Nodes->numNodes;  
   }  
   if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;  
   switch(elementtype)  
   {  
   case FINLEY_ELEMENTS:  
     elements=mesh_p->Elements;  
     break;  
   case FINLEY_FACE_ELEMENTS:  
     elements=mesh_p->FaceElements;  
     break;  
   case FINLEY_POINTS:  
     elements=mesh_p->Points;  
     break;  
   case FINLEY_CONTACT_ELEMENTS_1:  
     elements=mesh_p->ContactElements;  
     break;  
   }  
   if (elements==NULL)  
   {  
     Finley_setError(SYSTEM_ERROR,"saveVTK: undefined element file");  
     fclose(fileHandle_p);  
     TMPMEMFREE(isCellCentered);  
     return;  
   }  
   /* map finley element type to VTK element type */  
   numCells = elements->numElements;  
   if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)  
   {  
     TypeId = elements->LinearReferenceElement->Type->TypeId;  
   }  
   else  
   {  
     TypeId = elements->ReferenceElement->Type->TypeId;  
   }  
   
   switch(TypeId)  
   {  
   case Point1:  
   case Line2Face:  
   case Line3Face:  
   case Point1_Contact:  
   case Line2Face_Contact:  
   case Line3Face_Contact:  
     cellType = VTK_VERTEX;  
     numVTKNodesPerElement = 1;  
     strcpy(elemTypeStr, "VTK_VERTEX");  
     break;  
   
   case Line2:  
   case Tri3Face:  
   case Rec4Face:  
   case Line2_Contact:  
   case Tri3_Contact:  
   case Tri3Face_Contact:  
   case Rec4Face_Contact:  
     cellType = VTK_LINE;  
     numVTKNodesPerElement = 2;  
     strcpy(elemTypeStr, "VTK_LINE");  
     break;  
   
   case Tri3:  
   case Tet4Face:  
   case Tet4Face_Contact:  
     cellType = VTK_TRIANGLE;  
     numVTKNodesPerElement = 3;  
     strcpy(elemTypeStr, "VTK_TRIANGLE");  
     break;  
   
   case Rec4:  
   case Hex8Face:  
   case Rec4_Contact:  
   case Hex8Face_Contact:  
     cellType = VTK_QUAD;  
     numVTKNodesPerElement = 4;  
     strcpy(elemTypeStr, "VTK_QUAD");  
     break;  
   
   case Tet4:  
     cellType = VTK_TETRA;  
     numVTKNodesPerElement = 4;  
     strcpy(elemTypeStr, "VTK_TETRA");  
     break;  
   
   case Hex8:  
     cellType = VTK_HEXAHEDRON;  
     numVTKNodesPerElement = 8;  
     strcpy(elemTypeStr, "VTK_HEXAHEDRON");  
     break;  
   
   case Line3:  
   case Tri6Face:  
   case Rec8Face:  
   case Line3_Contact:  
   case Tri6Face_Contact:  
   case Rec8Face_Contact:  
     cellType = VTK_QUADRATIC_EDGE;  
     numVTKNodesPerElement = 3;  
     strcpy(elemTypeStr, "VTK_QUADRATIC_EDGE");  
     break;  
   
   case Tri6:  
   case Tet10Face:  
   case Tri6_Contact:  
   case Tet10Face_Contact:  
     cellType = VTK_QUADRATIC_TRIANGLE;  
     numVTKNodesPerElement = 6;  
     strcpy(elemTypeStr, "VTK_QUADRATIC_TRIANGLE");  
     break;  
   
   case Rec8:  
   case Hex20Face:  
   case Rec8_Contact:  
   case Hex20Face_Contact:  
     cellType = VTK_QUADRATIC_QUAD;  
     numVTKNodesPerElement = 8;  
     strcpy(elemTypeStr, "VTK_QUADRATIC_QUAD");  
     break;  
   
   case Tet10:  
     cellType = VTK_QUADRATIC_TETRA;  
     numVTKNodesPerElement = 10;  
     strcpy(elemTypeStr, "VTK_QUADRATIC_TETRA");  
     break;  
   
   case Hex20:  
     cellType = VTK_QUADRATIC_HEXAHEDRON;  
     numVTKNodesPerElement = 20;  
     strcpy(elemTypeStr, "VTK_QUADRATIC_HEXAHEDRON");  
     break;  
   
   default:  
     sprintf(error_msg, "saveVTK: Element type %s is not supported by VTK",elements->ReferenceElement->Type->Name);  
     Finley_setError(VALUE_ERROR,error_msg);  
     fclose(fileHandle_p);  
     TMPMEMFREE(isCellCentered);  
     return;  
   }  
   /* xml header */  
   fprintf(fileHandle_p, "<?xml version=\"1.0\"?>\n");  
   fprintf(fileHandle_p, "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\">\n");  
   
   /* finley uses an unstructured mesh, so UnstructuredGrid *should* work */  
   fprintf(fileHandle_p, "<UnstructuredGrid>\n");  
   
   /* is there only one "piece" to the data?? */  
   fprintf(fileHandle_p, "<Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n",numPoints, numCells);  
   /* now for the points; equivalent to positions section in saveDX() */  
   /* "The points element explicitly defines coordinates for each point  
   * individually.  It contains one DataArray element describing an array  
   * with three components per value, each specifying the coordinates of one  
   * point" - from Vtk User's Guide  
   */  
   fprintf(fileHandle_p, "<Points>\n");  
   /*  
   * the reason for this if statement is explained in the long comment below  
   */  
   nDim = mesh_p->Nodes->numDim;  
   fprintf(fileHandle_p, "<DataArray NumberOfComponents=\"%d\" type=\"Float64\" format=\"ascii\">\n",MAX(3,nDim));  
   /* vtk/mayavi doesn't like 2D data, it likes 3D data with a degenerate  
   * third dimension to handle 2D data (like a sheet of paper).  So, if  
   * nDim is 2, we have to append zeros to the array to get this third  
   * dimension, and keep the visualisers happy.  
   * Indeed, if nDim is less than 3, must pad all empty dimensions, so  
   * that the total number of dims is 3.  
   */  
   if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)  
   {  
     for (i = 0; i < mesh_p->Nodes->numNodes; i++)  
     {  
       if (mesh_p->Nodes->toReduced[i]>=0)  
       {  
         for (j = 0; j < nDim; j++) fprintf(fileHandle_p, " %e",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);  
         for (k=0; k<3-nDim; k++) fprintf(fileHandle_p, " %e",0.);  
         fprintf(fileHandle_p, "\n");  
       }  
     }  
   }  
   else  
   {  
     for (i = 0; i < mesh_p->Nodes->numNodes; i++)  
     {  
   
       for (j = 0; j < nDim; j++) fprintf(fileHandle_p, " %e",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);  
       for (k=0; k<3-nDim; k++) fprintf(fileHandle_p, " %e",0.);  
       fprintf(fileHandle_p, "\n");  
     }  
   
   }  
   fprintf(fileHandle_p, "</DataArray>\n");  
   fprintf(fileHandle_p, "</Points>\n");  
   
   /* write out the DataArray element for the connectivity */  
   
   NN = elements->ReferenceElement->Type->numNodes;  
   fprintf(fileHandle_p, "<Cells>\n");  
   fprintf(fileHandle_p, "<DataArray Name=\"connectivity\" type=\"Int32\" format=\"ascii\">\n");  
   
   if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)  
   {  
     for (i = 0; i < numCells; i++)  
     {  
       for (j = 0; j < numVTKNodesPerElement; j++)  
         fprintf(fileHandle_p,"%d ",mesh_p->Nodes->toReduced[elements->Nodes[INDEX2(elements->ReferenceElement->Type->linearNodes[j], i, NN)]]);  
       fprintf(fileHandle_p, "\n");  
     }  
   }  
   else if (VTK_QUADRATIC_HEXAHEDRON==cellType)  
   {  
     for (i = 0; i < numCells; i++)  
     {  
       fprintf(fileHandle_p,"%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d\n",  
               elements->Nodes[INDEX2(0, i, NN)],  
               elements->Nodes[INDEX2(1, i, NN)],  
               elements->Nodes[INDEX2(2, i, NN)],  
               elements->Nodes[INDEX2(3, i, NN)],  
               elements->Nodes[INDEX2(4, i, NN)],  
               elements->Nodes[INDEX2(5, i, NN)],  
               elements->Nodes[INDEX2(6, i, NN)],  
               elements->Nodes[INDEX2(7, i, NN)],  
               elements->Nodes[INDEX2(8, i, NN)],  
               elements->Nodes[INDEX2(9, i, NN)],  
               elements->Nodes[INDEX2(10, i, NN)],  
               elements->Nodes[INDEX2(11, i, NN)],  
               elements->Nodes[INDEX2(16, i, NN)],  
               elements->Nodes[INDEX2(17, i, NN)],  
               elements->Nodes[INDEX2(18, i, NN)],  
               elements->Nodes[INDEX2(19, i, NN)],  
               elements->Nodes[INDEX2(12, i, NN)],  
               elements->Nodes[INDEX2(13, i, NN)],  
               elements->Nodes[INDEX2(14, i, NN)],  
               elements->Nodes[INDEX2(15, i, NN)]);  
     }  
   }  
   else if (numVTKNodesPerElement!=NN)  
   {  
     for (i = 0; i < numCells; i++)  
     {  
       for (j = 0; j < numVTKNodesPerElement; j++) fprintf(fileHandle_p,"%d ", elements->Nodes[INDEX2(elements->ReferenceElement->Type->geoNodes[j], i, NN)]);  
       fprintf(fileHandle_p, "\n");  
     }  
   }  
   else  
   {  
     for (i = 0; i < numCells; i++)  
     {  
       for (j = 0; j < numVTKNodesPerElement; j++) fprintf(fileHandle_p,"%d ", elements->Nodes[INDEX2(j, i, NN)]);  
       fprintf(fileHandle_p, "\n");  
     }  
   }  
   fprintf(fileHandle_p, "</DataArray>\n");  
   
   /* write out the DataArray element for the offsets */  
   fprintf(fileHandle_p, "<DataArray Name=\"offsets\" type=\"Int32\" format=\"ascii\">\n");  
   for (i=numVTKNodesPerElement; i<=numCells*numVTKNodesPerElement; i+=numVTKNodesPerElement) fprintf(fileHandle_p, "%d\n", i);  
   fprintf(fileHandle_p, "</DataArray>\n");  
   
   /* write out the DataArray element for the types */  
   fprintf(fileHandle_p, "<DataArray Name=\"types\" type=\"UInt8\" format=\"ascii\">\n");  
   for (i=0; i<numCells; i++) fprintf(fileHandle_p, "%d\n", cellType);  
   fprintf(fileHandle_p, "</DataArray>\n");  
   
   /* finish off the <Cells> element */  
   fprintf(fileHandle_p, "</Cells>\n");  
   
   /* cell data */  
   if (write_celldata)  
   {  
     /* mark the active data arrays */  
     bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;  
     fprintf(fileHandle_p, "<CellData");  
     for (i_data =0 ;i_data<num_data;++i_data)  
     {  
       if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data])  
       {  
         /* if the rank == 0:   --> scalar data  
         * if the rank == 1:   --> vector data  
         * if the rank == 2:   --> tensor data  
         */  
         switch(getDataPointRank(data_pp[i_data]))  
         {  
         case 0:  
           if (! set_scalar)  
           {  
             fprintf(fileHandle_p," Scalars=\"%s\"",names_p[i_data]);  
             set_scalar=TRUE;  
           }  
           break;  
         case 1:  
           if (! set_vector)  
           {  
             fprintf(fileHandle_p," Vectors=\"%s\"",names_p[i_data]);  
             set_vector=TRUE;  
           }  
           break;  
         case 2:  
           if (! set_tensor)  
           {  
             fprintf(fileHandle_p," Tensors=\"%s\"",names_p[i_data]);  
             set_tensor=TRUE;  
           }  
           break;  
         default:  
           sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);  
           Finley_setError(VALUE_ERROR,error_msg);  
           fclose(fileHandle_p);  
           TMPMEMFREE(isCellCentered);  
           return;  
         }  
       }  
     }  
     fprintf(fileHandle_p, ">\n");  
     /* write the arrays */  
     for (i_data =0 ;i_data<num_data;++i_data)  
     {  
       if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data])  
       {  
         if (Finley_Assemble_reducedIntegrationOrder(data_pp[i_data])) {  
            numPointsPerSample=elements->ReferenceElementReducedOrder->numQuadNodes;  
         } else {  
            numPointsPerSample=elements->ReferenceElement->numQuadNodes;  
         }  
         rank = getDataPointRank(data_pp[i_data]);  
         nComp = getDataPointSize(data_pp[i_data]);  
         nCompReqd=1;   /* the number of components required by vtk */  
         shape=0;  
         if (rank == 0)  
         {  
           nCompReqd = 1;  
         }  
         else if (rank == 1)  
         {  
           shape=getDataPointShape(data_pp[i_data], 0);  
           if  (shape>3)  
           {  
             Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");  
             fclose(fileHandle_p);  
             TMPMEMFREE(isCellCentered);  
             return;  
           }  
           nCompReqd = 3;  
         }  
         else  
         {  
           shape=getDataPointShape(data_pp[i_data], 0);  
           if  (shape>3 || shape != getDataPointShape(data_pp[i_data], 1))  
           {  
             Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");  
             fclose(fileHandle_p);  
             TMPMEMFREE(isCellCentered);  
             return;  
           }  
           nCompReqd = 9;  
         }  
         fprintf(fileHandle_p, "<DataArray Name=\"%s\" type=\"Float64\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);  
   
         for (i=0; i<numCells; i++)  
         {  
           values = getSampleData(data_pp[i_data], i);  
           /* averaging over the number of points in the sample */  
           for (k=0; k<MIN(nComp,NCOMP_MAX); k++)  
           {  
             if (isExpanded(data_pp[i_data])) {  
                rtmp = 0.;  
                for (j=0; j<numPointsPerSample; j++) rtmp += values[INDEX2(k,j,nComp)];  
                sampleAvg[k] = rtmp/numPointsPerSample;  
867              } else {              } else {
868                 sampleAvg[k] = values[k];                shape=getDataPointShape(data_pp[i_data], 0);
869              }                if  (shape>3 || shape != getDataPointShape(data_pp[i_data], 1)) {
870                    Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
871            }                }
872            /* if the number of required components is more than the number                nCompReqd = 9;
873            * of actual components, pad with zeros              }
874            */              if (Finley_noError()) {
875            /* probably only need to get shape of first element */                 sprintf(txt_buffer,tag_Float_DataArray,names_p[i_data], nCompReqd);
876            /* write the data different ways for scalar, vector and tensor */                 if ( mpi_size > 1) {
877            if (nCompReqd == 1)                   if ( my_mpi_rank == 0) {
878            {                      #ifdef PASO_MPI
879              fprintf(fileHandle_p, " %e", sampleAvg[0]);                         MPI_File_iwrite_shared(mpi_fileHandle_p,txt_buffer,strlen(txt_buffer),MPI_CHAR,&mpi_req);
880            }                         MPI_Wait(&mpi_req,&mpi_status);
881            else if (nCompReqd == 3)                      #endif
882            {                   }
883              /* write out the data */                 } else {
884              for (m=0; m<shape; m++) fprintf(fileHandle_p, " %e", sampleAvg[m]);                     fprintf(fileHandle_p,txt_buffer);
885              for (m=0; m<nCompReqd-shape; m++) fprintf(fileHandle_p, " %e", 0.);                 }
886                   for (i=0; i<mesh_p->Nodes->numNodes; i++) {
887                      k=globalNodeIndex[i];
888                      if ( (myFirstNode <= k) && (k < myLastNode) ) {
889                         values = getSampleData(data_pp[i_data], nodeMapping->map[i]);
890                         /* if the number of mpi_required components is more than the number
891                         * of actual components, pad with zeros
892                         */
893                         /* probably only need to get shape of first element */
894                         /* write the data different ways for scalar, vector and tensor */
895                         if (nCompReqd == 1) {
896                           sprintf(tmp_buffer,FLOAT_SCALAR_FORMAT,values[0]);
897                         } else if (nCompReqd == 3) {
898                           if (shape==1) {
899                            sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,values[0],0.,0.);
900                           } else if (shape==2) {
901                            sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,values[0],values[1],0.);
902                           } else if (shape==3) {
903                            sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,values[0],values[2],values[3]);
904                           }
905                         } else if (nCompReqd == 9) {
906                           if (shape==1) {
907                            sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,values[0],0.,0.,
908                                                                   0.,0.,0.,
909                                                                   0.,0.,0.);
910                           } else if (shape==2) {
911                            sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,values[0],values[1],0.,
912                                                                   values[2],values[3],0.,
913                                                                   0.,0.,0.);
914                           } else if (shape==3) {
915                            sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,values[0],values[1],values[2],
916                                                                   values[3],values[4],values[5],
917                                                                   values[6],values[7],values[8]);
918                           }
919                         }
920                         if ( mpi_size > 1) {
921                           __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use);
922                         } else {
923                           fprintf(fileHandle_p,tmp_buffer);
924                         }
925                      }
926                   }
927                   if ( mpi_size > 1) {
928                       #ifdef PASO_MPI
929                         MPI_File_write_ordered(mpi_fileHandle_p,txt_buffer,txt_buffer_in_use, MPI_CHAR, &mpi_status);
930                       #endif    
931                       if ( my_mpi_rank == 0) {
932                          #ifdef PASO_MPI
933                             MPI_File_iwrite_shared(mpi_fileHandle_p,tag_End_DataArray,strlen(tag_End_DataArray),MPI_CHAR,&mpi_req);
934                             MPI_Wait(&mpi_req,&mpi_status);
935                          #endif
936                       }
937                   } else {
938                      fprintf(fileHandle_p,tag_End_DataArray);
939                   }
940                }
941              }
942            }
943            if ( mpi_size > 1) {
944              if ( my_mpi_rank == 0) {
945                 #ifdef PASO_MPI
946                    MPI_File_iwrite_shared(mpi_fileHandle_p,tag_End_CellData,strlen(tag_End_PointData),MPI_CHAR,&mpi_req);
947                    MPI_Wait(&mpi_req,&mpi_status);
948                 #endif
949            }            }
           else if (nCompReqd == 9)  
           {  
             /* tensor data, so have a 3x3 matrix to output as a row  
             * of 9 data points */  
             count = 0;  
             for (m=0; m<shape; m++)  
             {  
               for (n=0; n<shape; n++)  
               {  
                 fprintf(fileHandle_p, " %e", sampleAvg[count]);  
                 count++;  
               }  
               for (n=0; n<3-shape; n++) fprintf(fileHandle_p, " %e", 0.);  
             }  
             for (m=0; m<3-shape; m++)  
               for (n=0; n<3; n++) fprintf(fileHandle_p, " %e", 0.);  
           }  
           fprintf(fileHandle_p, "\n");  
         }  
         fprintf(fileHandle_p, "</DataArray>\n");  
       }  
     }  
     fprintf(fileHandle_p, "</CellData>\n");  
   }  
   /* point data */  
   if (write_pointdata)  
   {  
     /* mark the active data arrays */  
     bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;  
     fprintf(fileHandle_p, "<PointData");  
     for (i_data =0 ;i_data<num_data;++i_data)  
     {  
       if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])  
       {  
         /* if the rank == 0:   --> scalar data  
         * if the rank == 1:   --> vector data  
         * if the rank == 2:   --> tensor data  
         */  
         switch(getDataPointRank(data_pp[i_data]))  
         {  
         case 0:  
           if (! set_scalar)  
           {  
             fprintf(fileHandle_p," Scalars=\"%s\"",names_p[i_data]);  
             set_scalar=TRUE;  
           }  
           break;  
         case 1:  
           if (! set_vector)  
           {  
             fprintf(fileHandle_p," Vectors=\"%s\"",names_p[i_data]);  
             set_vector=TRUE;  
           }  
           break;  
         case 2:  
           if (! set_tensor)  
           {  
             fprintf(fileHandle_p," Tensors=\"%s\"",names_p[i_data]);  
             set_tensor=TRUE;  
           }  
           break;  
         default:  
           sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);  
           Finley_setError(VALUE_ERROR,error_msg);  
           fclose(fileHandle_p);  
           TMPMEMFREE(isCellCentered);  
           return;  
         }  
       }  
     }  
     fprintf(fileHandle_p, ">\n");  
     /* write the arrays */  
     for (i_data =0 ;i_data<num_data;++i_data)  
     {  
       if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])  
       {  
         if (Finley_Assemble_reducedIntegrationOrder(data_pp[i_data])) {  
            numPointsPerSample=elements->ReferenceElementReducedOrder->numQuadNodes;  
950          } else {          } else {
951             numPointsPerSample=elements->ReferenceElement->numQuadNodes;              fprintf(fileHandle_p,tag_End_PointData);
         }  
         rank = getDataPointRank(data_pp[i_data]);  
         nComp = getDataPointSize(data_pp[i_data]);  
         nCompReqd=1;   /* the number of components required by vtk */  
         shape=0;  
         if (rank == 0)  
         {  
           nCompReqd = 1;  
952          }          }
         else if (rank == 1)  
         {  
           shape=getDataPointShape(data_pp[i_data], 0);  
           if  (shape>3)  
           {  
             Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");  
             fclose(fileHandle_p);  
             TMPMEMFREE(isCellCentered);  
             return;  
           }  
           nCompReqd = 3;  
         }  
         else  
         {  
           shape=getDataPointShape(data_pp[i_data], 0);  
           if  (shape>3 || shape != getDataPointShape(data_pp[i_data], 1))  
           {  
             Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");  
             fclose(fileHandle_p);  
             TMPMEMFREE(isCellCentered);  
             return;  
           }  
           nCompReqd = 9;  
         }  
         fprintf(fileHandle_p, "<DataArray Name=\"%s\" type=\"Float64\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);  
         /* write out the data */  
         /* if the number of required components is more than the number  
         * of actual components, pad with zeros  
         */  
         do_write=TRUE;  
         for (i=0; i<mesh_p->Nodes->numNodes; i++)  
         {  
           if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)  
           {  
             if (mesh_p->Nodes->toReduced[i]>=0)  
             {  
               switch(getFunctionSpaceType(data_pp[i_data]))  
               {  
               case FINLEY_DEGREES_OF_FREEDOM:  
                 values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);  
                 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;  
             }  
           }  
           else  
           {  
             do_write=TRUE;  
             switch(getFunctionSpaceType(data_pp[i_data]))  
             {  
             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;  
             }  
           }  
           /* write the data different ways for scalar, vector and tensor */  
           if (do_write)  
           {  
             if (nCompReqd == 1)  
             {  
               fprintf(fileHandle_p, " %e", values[0]);  
             }  
             else if (nCompReqd == 3)  
             {  
               for (m=0; m<shape; m++) fprintf(fileHandle_p, " %e", values[m]);  
               for (m=0; m<nCompReqd-shape; m++) fprintf(fileHandle_p, " %e", 0.);  
             }  
             else if (nCompReqd == 9)  
             {  
               /* tensor data, so have a 3x3 matrix to output as a row  
               * of 9 data points */  
               count = 0;  
               for (m=0; m<shape; m++)  
               {  
                 for (n=0; n<shape; n++)  
                 {  
                   fprintf(fileHandle_p, " %e", values[count]);  
                   count++;  
                 }  
                 for (n=0; n<3-shape; n++) fprintf(fileHandle_p, " %e", 0.);  
               }  
               for (m=0; m<3-shape; m++)  
                 for (n=0; n<3; n++) fprintf(fileHandle_p, " %e", 0.);  
             }  
             fprintf(fileHandle_p, "\n");  
           }  
         }  
         fprintf(fileHandle_p, "</DataArray>\n");  
       }  
     }  
     fprintf(fileHandle_p, "</PointData>\n");  
953    }    }
954    /* finish off the piece */    if (Finley_noError()) {
955    fprintf(fileHandle_p, "</Piece>\n");       if ( mpi_size > 1) {
956           if ( my_mpi_rank == 0) {
957    fprintf(fileHandle_p, "</UnstructuredGrid>\n");            #ifdef PASO_MPI
958    /* write the xml footer */               MPI_File_iwrite_shared(mpi_fileHandle_p,footer,strlen(footer),MPI_CHAR,&mpi_req);
959    fprintf(fileHandle_p, "</VTKFile>\n");               MPI_Wait(&mpi_req,&mpi_status);
960    /* close the file */               #ifdef MPIO_HINTS
961    fclose(fileHandle_p);                 MPI_Info_free(&mpi_info);
962                   #undef MPIO_HINTS
963                 #endif
964                 MPI_File_close(&mpi_fileHandle_p);
965              #endif
966            }
967         } else {
968             fprintf(fileHandle_p,footer);
969             fclose(fileHandle_p);
970         }
971      }
972    TMPMEMFREE(isCellCentered);    TMPMEMFREE(isCellCentered);
973      TMPMEMFREE(txt_buffer);
974    return;    return;
975  }  }
 #endif  
   

Legend:
Removed from v.1222  
changed lines
  Added in v.1223

  ViewVC Help
Powered by ViewVC 1.1.26