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

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

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

revision 903 by gross, Fri Nov 17 01:59:49 2006 UTC revision 2094 by caltinay, Tue Nov 25 05:45:41 2008 UTC
# Line 1  Line 1 
 /*  
 ************************************************************  
 *          Copyright 2006 by ACcESS MNRF                   *  
 *                                                          *  
 *              http://www.access.edu.au                    *  
 *       Primary Business: Queensland, Australia            *  
 *  Licensed under the Open Software License version 3.0    *  
 *     http://www.opensource.org/licenses/osl-3.0.php       *  
 *                                                          *  
 ************************************************************  
 */  
1    
2    /*******************************************************
3    *
4    * Copyright (c) 2003-2008 by University of Queensland
5    * Earth Systems Science Computational Center (ESSCC)
6    * http://www.uq.edu.au/esscc
7    *
8    * Primary Business: Queensland, Australia
9    * Licensed under the Open Software License version 3.0
10    * http://www.opensource.org/licenses/osl-3.0.php
11    *
12    *******************************************************/
13    
 /**************************************************************/  
   
 /*   writes data and mesh in a vtk file */  
14    
15  /**************************************************************/  /**************************************************************/
16    
17  /*   Author: Paul Cochrane, cochrane@esscc.uq.edu.au */  /*   writes data and mesh in a vtk file */
18  /*   MPI-IO version: Derick Hawcroft, d.hawcroft@uq.edu.au         */  /*   nodal data needs to be given on FINLEY_NODES or FINLEY_REDUCED_NODES */
   
 /*   Version: $Id$ */  
19    
20  /**************************************************************/  /**************************************************************/
21    
22    
23  #include "Mesh.h"  #include "Mesh.h"
24    #include "Assemble.h"
25  #include "vtkCellType.h"  /* copied from vtk source directory !!! */  #include "vtkCellType.h"  /* copied from vtk source directory !!! */
26    #include "paso/PasoUtil.h"
27    
28  /*  #define LEN_PRINTED_INT_FORMAT (9+1)
29   MPI version notes:  #define INT_FORMAT "%d "
30    #define INT_NEWLINE_FORMAT "%d\n"
31   ******************************************************************************  #define FLOAT_SCALAR_FORMAT "%12.6e\n"
32   ***                                                                       ****  #define FLOAT_VECTOR_FORMAT "%12.6e %12.6e %12.6e\n"
33   *** WARNING: Won't work for meshes with peridodic boundary conditions yet ****  #define FLOAT_TENSOR_FORMAT "%12.6e %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e\n"
34   ***                                                                       ****    #define LEN_PRINTED_FLOAT_SCALAR_FORMAT (12+1)
35   ******************************************************************************  #define LEN_PRINTED_FLOAT_VECTOR_FORMAT (3*(12+1)+1)
36    #define LEN_PRINTED_FLOAT_TENSOR_FORMAT (9*(12+1)+1)
37   In this version, the rank==0 process writes *all* opening and closing  #define NEWLINE "\n"
38   XML tags.  #define LEN_TMP_BUFFER LEN_PRINTED_FLOAT_TENSOR_FORMAT+(MAX_numNodes*LEN_PRINTED_INT_FORMAT+1)+2
39   Individual process data is copied to a buffer before being written  #define NCOMP_MAX 9
40   out. The  routines are collectively called and will be called in the natural  #define __STRCAT(dest,chunk,dest_in_use)  \
41   ordering i.e 0 to maxProcs-1.  {                  \
42      strcpy(&dest[dest_in_use], chunk); \
43   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); \  
44  }  }
45    #define INSIDE_1D(_X_,_C_,_R_) ( ABS((_X_)-(_C_)) <= (_R_) )
46  void Finley_Mesh_saveVTK_MPIO(const char * filename_p, Finley_Mesh *mesh_p, const dim_t num_data,char* *names_p, escriptDataC* *data_pp)  #define INSIDE_2D(_X_,_Y_,_CX_,_CY_,_R_) ( INSIDE_1D(_X_,_CX_,_R_) &&  INSIDE_1D(_Y_,_CY_,_R_))
47    #define INSIDE_3D(_X_,_Y_,_Z_,_CX_,_CY_,_CZ_,_R_) ( INSIDE_1D(_X_,_CX_,_R_) &&  INSIDE_1D(_Y_,_CY_,_R_) && INSIDE_1D(_Z_,_CZ_,_R_) )
48    
49    void Finley_Mesh_saveVTK(const char * filename_p,
50                             Finley_Mesh *mesh_p,
51                             const dim_t num_data,
52                             char* *names_p,
53                             escriptDataC* *data_pp)
54  {  {
55    int    numPoints,  #ifdef USE_VTK
56    numCells = -1,    char error_msg[LenErrorMsg_MAX], *txt_buffer=NULL, tmp_buffer[LEN_TMP_BUFFER];
57               myRank,comm,gsize,    double sampleAvg[NCOMP_MAX], *values, *QuadNodes=NULL;
58               numLocal,    size_t txt_buffer_in_use;
59               nDim,    dim_t len_txt_buffer,  max_len_names;
60               shape;    FILE * fileHandle_p = NULL;
61    size_t __n;    int mpi_size, i, j, l, cellType=0;
62    int i,j,k,m,n,count;    dim_t i_data, hits, hits_old;
63    int numGlobalCells = 0;    dim_t nDim, globalNumPoints=0, numCells=0, globalNumCells=0, numVTKNodesPerElement=0;
64    index_t  *nodesGlobal=NULL;   // used to get the connectivity  right for VTK    dim_t myNumPoints=0, numPointsPerSample, rank, nComp, nCompReqd;
65      dim_t shape, NN=0, numCellFactor=0, myNumCells=0, max_name_len;
66    /* variables associatted with write_celldata/pointdata */    bool_t *isCellCentered=NULL, write_celldata=FALSE, write_pointdata=FALSE, reduced_elements=FALSE;
67    int numPointsPerSample,    bool_t set_scalar=FALSE, set_vector=FALSE, set_tensor=FALSE;
68    nComp,    index_t myFirstNode=0, myLastNode=0, *globalNodeIndex=NULL;
69    nCompReqd;    index_t k, *node_index, myFirstCell=0;
70    double* values, rtmp;    #ifdef PASO_MPI
71      int ierr;
72    // Local element info (for debugging)    /* int amode = MPI_MODE_CREATE | MPI_MODE_WRONLY |  MPI_MODE_SEQUENTIAL;  */
73    size_t numLocalCells,    const int amode = MPI_MODE_CREATE | MPI_MODE_WRONLY | MPI_MODE_UNIQUE_OPEN;
74    numInternalCells,    MPI_File mpi_fileHandle_p;
75    numBoundaryCells;    MPI_Status mpi_status;
76      MPI_Request mpi_req;
77    int rank;    MPI_Info mpi_info=MPI_INFO_NULL;
78      #endif
79    int amode = MPI_MODE_CREATE | MPI_MODE_WRONLY |  MPI_MODE_SEQUENTIAL;    Paso_MPI_rank my_mpi_rank;
80      int nodetype=FINLEY_NODES;
   comm   = mesh_p->Nodes->MPIInfo->comm;  
   myRank = mesh_p->Nodes->MPIInfo->rank;  
   gsize  = mesh_p->Nodes->MPIInfo->size;  
   
   MPI_File fh;  
   MPI_Status status;  
   MPI_Request req;  
   MPI_Info infoHints;  
   
   char error_msg[LenErrorMsg_MAX];  
   
   int i_data;  
   
   int nodetype=FINLEY_DEGREES_OF_FREEDOM;  
81    int elementtype=FINLEY_UNKNOWN;    int elementtype=FINLEY_UNKNOWN;
82    bool_t isCellCentered[num_data],write_celldata=FALSE,write_pointdata=FALSE;    Finley_NodeMapping *nodeMapping=NULL;
   
   ElementTypeId TypeId;  
   
   int numVTKNodesPerElement;  
   int cellType;  
   char elemTypeStr[32];  
   
83    Finley_ElementFile* elements=NULL;    Finley_ElementFile* elements=NULL;
84      ElementTypeId TypeId=NoType;
85      
86    
87      /****************************************/
88      /*                                      */
89      /*       tags in the vtk file           */
90    
91      const char* tags_header="<?xml version=\"1.0\"?>\n" \
92                        "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\">\n" \
93                        "<UnstructuredGrid>\n" \
94                        "<Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n" \
95                        "<Points>\n" \
96                        "<DataArray NumberOfComponents=\"%d\" type=\"Float64\" format=\"ascii\">\n";
97      char* tag_End_DataArray = "</DataArray>\n";
98      char* tag_End_PointData = "</PointData>\n";
99      char* tag_End_CellData =  "</CellData>\n";
100      char* footer = "</Piece>\n</UnstructuredGrid>\n</VTKFile>\n";
101      char* tags_End_Points_and_Start_Conn = "</DataArray>\n</Points>\n<Cells>\n<DataArray Name=\"connectivity\" type=\"Int32\" format=\"ascii\">\n" ;
102      char* tags_End_Conn_and_Start_Offset = "</DataArray>\n<DataArray Name=\"offsets\" type=\"Int32\" format=\"ascii\">\n";
103      char* tags_End_Offset_and_Start_Type = "</DataArray>\n<DataArray Name=\"types\" type=\"UInt8\" format=\"ascii\">\n";
104      const char* tag_Float_DataArray="<DataArray Name=\"%s\" type=\"Float64\" NumberOfComponents=\"%d\" format=\"ascii\">\n";
105      char* tags_End_Type_And_Cells = "</DataArray>\n</Cells>\n";
106    
107      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 };
108    
109    // Local node info    /* if there is no mesh we just return */
110    int numInternalNodes,    if (mesh_p==NULL) return;
   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:  
         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:  
         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:  
         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:  
         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 ;  
   }  
111    
112    numCells =  elements->numElements;    my_mpi_rank = mesh_p->Nodes->MPIInfo->rank;
113    numGlobalCells = elements->elementDistribution->vtxdist[gsize];    mpi_size  = mesh_p->Nodes->MPIInfo->size;
114    numLocalCells    = elements->elementDistribution->numLocal;    nDim = mesh_p->Nodes->numDim;
   numInternalCells = elements->elementDistribution->numInternal;  
   numBoundaryCells = elements->elementDistribution->numBoundary;  
   
   if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)  
   {  
     TypeId = elements->LinearReferenceElement->Type->TypeId;  
   }  
   else  
   {  
     TypeId = elements->ReferenceElement->Type->TypeId;  
   }  
115    
116    switch(TypeId)    if (! ( (nDim ==2) || (nDim == 3) ) ) {
117    {          Finley_setError(IO_ERROR, "saveVTK: spatial dimension 2 or 3 is supported only.");
118    case Point1:          return;  
   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;  
119    }    }
120      /*************************************************************************************/
121    
122    /* Write XML Header */    /* open the file and check handle */
   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));  
   
123    
124      MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);    if (mpi_size > 1) {
125      MPI_Wait(&req,&status);          #ifdef PASO_MPI
126              /* Collective Call */
127              #ifdef MPIO_HINTS
128                MPI_Info_create(&mpi_info);
129                /*  MPI_Info_set(mpi_info, "striping_unit",        "424288"); */
130                /*  MPI_Info_set(mpi_info, "striping_factor",      "16"); */
131                /*  MPI_Info_set(mpi_info, "collective_buffering", "true"); */
132                /*  MPI_Info_set(mpi_info, "cb_block_size",        "131072"); */
133                /*  MPI_Info_set(mpi_info, "cb_buffer_size",       "1048567"); */
134                /*  MPI_Info_set(mpi_info, "cb_nodes",             "8"); */
135                /*    MPI_Info_set(mpi_info, "access_style", "write_once, sequential"); */
136              
137                /*XFS only */
138                /*   MPI_Info_set(mpi_info, "direct_write",          "true"); */
139              #endif
140              if ( my_mpi_rank == 0) {
141                  if  (Paso_fileExists(filename_p)) remove(filename_p);
142              }
143              ierr=MPI_File_open(mesh_p->Nodes->MPIInfo->comm, (char*)filename_p, amode,mpi_info, &mpi_fileHandle_p);
144              if (ierr != MPI_SUCCESS) {
145              perror(filename_p);
146                  sprintf(error_msg, "saveVTK: File %s could not be opened for writing in parallel.", filename_p);
147                  Finley_setError(IO_ERROR,error_msg);
148              } else {
149                 MPI_File_set_view(mpi_fileHandle_p,MPI_DISPLACEMENT_CURRENT,MPI_CHAR, MPI_CHAR, "native" , mpi_info);
150              }
151            #endif
152      } else {
153            fileHandle_p = fopen(filename_p, "w");
154            if (fileHandle_p==NULL) {
155               sprintf(error_msg, "saveVTK: File %s could not be opened for writing.", filename_p);
156               Finley_setError(IO_ERROR,error_msg);
157             }
158    }    }
159      if (! Paso_MPIInfo_noError(mesh_p->Nodes->MPIInfo) ) return;
160      /*************************************************************************************/
161    
162    MPIO_DEBUG(" Writing Coordinate Points... ")    /* find the mesh type to be written */
   
   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;              \  
 }  
     
   // Loop over all nodes      
   for (i = 0; i < mesh_p->Nodes->numNodes; i++)  
   {  
     // This is the bit that will break for periodic BCs because it assumes that there is a one to one  
     // correspondance between nodes and Degrees of freedom  
     //TODO: handle periodic BC's  
     DOFNodes[mesh_p->Nodes->degreeOfFreedom[i]] = i;  
   
     // Is this node local to the domain ?  
     if( mesh_p->Nodes->degreeOfFreedom[i] < localDOF )  
     {  
       for (j = 0; j < nDim; j++)  
       {  
         sprintf(tmpbuf,"%e ", mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)] );  
         __STRCAT(largebuf,tmpbuf,tsz)  
       }  
       for (k=0; k<3-nDim; k++)  
       {  
         __STRCAT(largebuf,zero,tsz)  
       }  
       __STRCAT(largebuf,newline,tsz)  
       nodeCache.values[numNodesOutput++]=i;  
     }  
   }  
163    
164    nodeCache.size=numNodesOutput;    isCellCentered=TMPMEMALLOC(num_data,bool_t);
165      max_len_names=0;
166      if (!Finley_checkPtr(isCellCentered)) {
167         reduced_elements=FALSE;
168         nodetype=FINLEY_UNKNOWN;
169         elementtype=FINLEY_UNKNOWN;
170         for (i_data=0;i_data<num_data;++i_data) {
171           if (! isEmpty(data_pp[i_data])) {
172             switch(getFunctionSpaceType(data_pp[i_data]) ) {
173             case FINLEY_NODES:
174               nodetype = (nodetype == FINLEY_REDUCED_NODES) ? FINLEY_REDUCED_NODES : FINLEY_NODES;
175               isCellCentered[i_data]=FALSE;
176               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
177                 elementtype=FINLEY_ELEMENTS;
178               } else {
179                 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
180               }
181               break;
182             case FINLEY_REDUCED_NODES:
183               nodetype = FINLEY_REDUCED_NODES;
184               isCellCentered[i_data]=FALSE;
185               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
186                 elementtype=FINLEY_ELEMENTS;
187               } else {
188                 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
189               }
190               break;
191             case FINLEY_REDUCED_ELEMENTS:
192                reduced_elements=TRUE;
193             case FINLEY_ELEMENTS:
194               isCellCentered[i_data]=TRUE;
195               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
196                 elementtype=FINLEY_ELEMENTS;
197               } else {
198                 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
199               }
200               break;
201             case FINLEY_REDUCED_FACE_ELEMENTS:
202                reduced_elements=TRUE;
203             case FINLEY_FACE_ELEMENTS:
204               isCellCentered[i_data]=TRUE;
205               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_FACE_ELEMENTS) {
206                 elementtype=FINLEY_FACE_ELEMENTS;
207               } else {
208                 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
209               }
210               break;
211             case FINLEY_POINTS:
212               isCellCentered[i_data]=TRUE;
213               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_POINTS) {
214                 elementtype=FINLEY_POINTS;
215               } else {
216                 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
217               }
218               break;
219             case FINLEY_REDUCED_CONTACT_ELEMENTS_1:
220                reduced_elements=TRUE;
221             case FINLEY_CONTACT_ELEMENTS_1:
222               isCellCentered[i_data]=TRUE;
223               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {
224                 elementtype=FINLEY_CONTACT_ELEMENTS_1;
225               } else {
226                 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
227               }
228               break;
229             case FINLEY_REDUCED_CONTACT_ELEMENTS_2:
230                reduced_elements=TRUE;
231             case FINLEY_CONTACT_ELEMENTS_2:
232               isCellCentered[i_data]=TRUE;
233               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {
234                 elementtype=FINLEY_CONTACT_ELEMENTS_1;
235               } else {
236                 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
237               }
238               break;
239             default:
240               sprintf(error_msg,"saveVTK: unknown function space type %d",getFunctionSpaceType(data_pp[i_data]));
241               Finley_setError(TYPE_ERROR,error_msg);
242             }
243             if (isCellCentered[i_data]) {
244               write_celldata=TRUE;
245             } else {
246               write_pointdata=TRUE;
247             }
248             max_len_names =MAX(max_len_names,(dim_t)strlen(names_p[i_data]));
249           }
250         }
251         nodetype = (nodetype == FINLEY_UNKNOWN) ? FINLEY_NODES : nodetype;
252      }
253      if (Finley_noError()) {
254    
255         /***************************************/
256    
257         /* select number of points and the mesh component */
258    
259         if (nodetype == FINLEY_REDUCED_NODES) {
260            myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
261            myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
262            globalNumPoints = Finley_NodeFile_getGlobalNumReducedNodes(mesh_p->Nodes);
263            globalNodeIndex= Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
264         } else {
265            myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
266            myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
267            globalNumPoints = Finley_NodeFile_getGlobalNumNodes(mesh_p->Nodes);
268            globalNodeIndex= Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
269         }
270         myNumPoints = myLastNode - myFirstNode;
271         if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;
272         switch(elementtype) {
273           case FINLEY_ELEMENTS:
274              elements=mesh_p->Elements;
275              break;
276            case FINLEY_FACE_ELEMENTS:
277              elements=mesh_p->FaceElements;
278              break;
279            case FINLEY_POINTS:
280              elements=mesh_p->Points;
281              break;
282            case FINLEY_CONTACT_ELEMENTS_1:
283              elements=mesh_p->ContactElements;
284              break;
285         }
286         if (elements==NULL) {
287           Finley_setError(SYSTEM_ERROR,"saveVTK: undefined element file");
288         } else {
289           /* map finley element type to VTK element type */
290           numCells = elements->numElements;
291           globalNumCells = Finley_ElementFile_getGlobalNumElements(elements);
292           myNumCells= Finley_ElementFile_getMyNumElements(elements);
293           myFirstCell= Finley_ElementFile_getFirstElement(elements);
294           NN = elements->numNodes;
295           if (nodetype==FINLEY_REDUCED_NODES) {
296              TypeId = elements->LinearReferenceElement->Type->TypeId;
297              if (reduced_elements) {
298                  QuadNodes=elements->LinearReferenceElementReducedOrder->QuadNodes;
299              } else {
300                  QuadNodes=elements->LinearReferenceElement->QuadNodes;
301              }
302           } else {
303              TypeId = elements->ReferenceElement->Type->TypeId;
304              if (reduced_elements) {
305                  QuadNodes=elements->ReferenceElementReducedOrder->QuadNodes;
306              } else {
307                  QuadNodes=elements->ReferenceElement->QuadNodes;
308              }
309           }
310           switch(TypeId) {
311            case Point1:
312            case Line2Face:
313            case Line3Face:
314            case Point1_Contact:
315            case Line2Face_Contact:
316            case Line3Face_Contact:
317              numCellFactor=1;
318              cellType = VTK_VERTEX;
319              numVTKNodesPerElement = 1;
320              break;
321          
322            case Line2:
323            case Tri3Face:
324            case Rec4Face:
325            case Line2_Contact:
326            case Tri3_Contact:
327            case Tri3Face_Contact:
328            case Rec4Face_Contact:
329              numCellFactor=1;
330              cellType = VTK_LINE;
331              numVTKNodesPerElement = 2;
332              break;
333          
334            case Tri3:
335            case Tet4Face:
336            case Tet4Face_Contact:
337              numCellFactor=1;
338              cellType = VTK_TRIANGLE;
339              numVTKNodesPerElement = 3;
340              break;
341          
342            case Rec4:
343            case Hex8Face:
344            case Rec4_Contact:
345            case Hex8Face_Contact:
346              numCellFactor=1;
347              cellType = VTK_QUAD;
348              numVTKNodesPerElement = 4;
349              break;
350    
351    largebuf[tsz] = '\0';          case Rec9:
352    MPI_File_write_ordered(fh, largebuf,tsz, MPI_CHAR, &status);            numCellFactor=4;
353    MEMFREE(largebuf);            cellType = VTK_QUAD;
354              numVTKNodesPerElement = 4;
355    nodesGlobal = MEMALLOC(mesh_p->Nodes->numNodes ,index_t);            break;
356          
357    // form distribution info on who output which nodes          case Tet4:
358    vtxdist = MEMALLOC( gsize+1, index_t );            numCellFactor=1;
359    vtxdist[0]=0;            cellType = VTK_TETRA;
360    MPI_Allgather(&numNodesOutput,1,MPI_INT,vtxdist+1,1,MPI_INT,comm);            numVTKNodesPerElement = 4;
361    for( i=0; i<gsize; i++ )            break;
362      vtxdist[i+1]+=vtxdist[i];        
363            case Hex8:
364    // will not work for periodic boundary conditions            numCellFactor=1;
365    // calculate the local nodes file positions            cellType = VTK_HEXAHEDRON;
366    pos = 0;            numVTKNodesPerElement = 8;
367    for( i=0; i<mesh_p->Nodes->numNodes; i++ )            break;
368    {        
369      if( mesh_p->Nodes->degreeOfFreedom[i]< localDOF )          case Line3:
370      {          case Tri6Face:
371        nodesGlobal[i] = vtxdist[myRank] + pos++;          case Rec8Face:
372      }          case Line3_Contact:
373      else          case Tri6Face_Contact:
374        nodesGlobal[i] = -1;          case Rec8Face_Contact:
375    }            numCellFactor=1;
376              cellType = VTK_QUADRATIC_EDGE;
377              numVTKNodesPerElement = 3;
378              break;
379          
380            case Tri6:
381            case Tet10Face:
382            case Tri6_Contact:
383            case Tet10Face_Contact:
384              numCellFactor=1;
385              cellType = VTK_QUADRATIC_TRIANGLE;
386              numVTKNodesPerElement = 6;
387              break;
388          
389            case Rec8:
390            case Hex20Face:
391            case Rec8_Contact:
392            case Hex20Face_Contact:
393              numCellFactor=1;
394              cellType = VTK_QUADRATIC_QUAD;
395              numVTKNodesPerElement = 8;
396              break;
397          
398            case Tet10:
399              numCellFactor=1;
400              cellType = VTK_QUADRATIC_TETRA;
401              numVTKNodesPerElement = 10;
402              break;
403          
404            case Hex20:
405              numCellFactor=1;
406              cellType = VTK_QUADRATIC_HEXAHEDRON;
407              numVTKNodesPerElement = 20;
408              break;
409    
410    // communicate the local Nodes file position to the interested parties          case Hex27:
411    // send local info            numCellFactor=8;
412    forwardBuffer = MEMALLOC( mesh_p->Nodes->numNodes, index_t );            cellType = VTK_HEXAHEDRON;
413    for( n=0; n < dist->numNeighbours; n++ )            numVTKNodesPerElement = 8;
414    {            break;
415      if(  dist->edges[n]->numForward)        
416      {          default:
417        for( i=0; i < dist->edges[n]->numForward; i++ )            sprintf(error_msg, "saveVTK: Element type %s is not supported by VTK",elements->ReferenceElement->Type->Name);
418          forwardBuffer[i] = nodesGlobal[DOFNodes[dist->edges[n]->indexForward[i] ]];            Finley_setError(VALUE_ERROR,error_msg);
419        Paso_CommBuffer_pack( mesh_p->Nodes->CommBuffer, dist->neighbours[n], NULL, forwardBuffer, sizeof(index_t), 0 );          }
420        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];  
     }  
421    }    }
422        /***************************************/
423    
424      /***************************************/
425      /*                                     */
426      /*   allocate text buffer              */
427      /*                                     */
428      max_name_len=0;
429      for (i_data =0 ;i_data<num_data;++i_data) max_name_len=MAX(max_name_len,(dim_t)strlen(names_p[i_data]));
430      len_txt_buffer= strlen(tags_header) + 3 * LEN_PRINTED_INT_FORMAT + (30+3*max_name_len); /* header */
431      if (mpi_size > 1) len_txt_buffer=MAX(len_txt_buffer, myNumPoints * LEN_TMP_BUFFER);
432      if (mpi_size > 1) len_txt_buffer=MAX(len_txt_buffer, numCellFactor*myNumCells*(LEN_PRINTED_INT_FORMAT*numVTKNodesPerElement+1));
433      len_txt_buffer=MAX(len_txt_buffer,200+3*max_len_names);
434      len_txt_buffer=MAX(len_txt_buffer, (dim_t)strlen(tag_Float_DataArray) + LEN_PRINTED_INT_FORMAT + max_len_names);
435      if (mpi_size > 1) len_txt_buffer=MAX(len_txt_buffer, numCellFactor*myNumCells*LEN_PRINTED_FLOAT_TENSOR_FORMAT);
436      if (mpi_size > 1) len_txt_buffer=MAX(len_txt_buffer, myNumPoints*LEN_PRINTED_FLOAT_TENSOR_FORMAT);
437      txt_buffer=TMPMEMALLOC(len_txt_buffer+1,char);
438      Finley_checkPtr(txt_buffer);
439        
440    MEMFREE(vtxdist);    if (Finley_noError()) {
   MEMFREE(DOFNodes);  
   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);  
   }  
   MPIO_DEBUG(" Done Writing Coordinate Points ")  
   
   /* BEGIN CONNECTIVITY */  
441    
442    int NN = elements->ReferenceElement->Type->numNodes; /* num Nodes holding ref-element */       /* select number of points and the mesh component */
443    
444    // Collective       sprintf(txt_buffer,tags_header,globalNumPoints,numCellFactor*globalNumCells,3);
   MPIO_DEBUG(" Writing Connectivity... ")  
445    
446    size_t sz = numLocalCells*6*numVTKNodesPerElement + numLocalCells;        if (mpi_size > 1) {
447    largebuf = MEMALLOC(sz,char);            if ( my_mpi_rank == 0) {
448    largebuf[0] = '\0';              #ifdef PASO_MPI
449    tsz=0;                MPI_File_iwrite_shared(mpi_fileHandle_p,txt_buffer,strlen(txt_buffer),MPI_CHAR,&mpi_req);
450    pos = 0;                MPI_Wait(&mpi_req,&mpi_status);
451    // numCells?              #endif
452    elementCache.values = MEMALLOC(numLocalCells,index_t);            }
453    if (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM)        } else {
454    {           fprintf(fileHandle_p, "%s", txt_buffer);
455      for (i = 0; i < numCells; i++)        }
456      {  
457          /* write the nodes */
458        if (elements->Id[i] >= elements->elementDistribution->vtxdist[i] &&  elements->Id[i] <= elements->elementDistribution->vtxdist[i+1] - 1 )        
459        {        if (mpi_size > 1) {
460          for (j = 0; j < numVTKNodesPerElement; j++)  
461          {           txt_buffer[0] = '\0';
462            sprintf(tmpbuf,"%d ",nodesGlobal[mesh_p->Nodes->toReduced[elements->Nodes[INDEX2(elements->ReferenceElement->Type->linearNodes[j], i, NN)]]]);           txt_buffer_in_use=0;
463            __STRCAT(largebuf,tmpbuf,tsz)           if (nDim==2) {
464          }              for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
465          __STRCAT(largebuf,newline,tsz)                 if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
466          elementCache.values[pos++]=i;                   sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,
467        }                                      mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
468      }                                      mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
469    }                                      0.);
470    else if (VTK_QUADRATIC_HEXAHEDRON==cellType)                   __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use);
471    {                 }
472      char tmpbuf2[20*20*2];              }      
473             } else {
474      for (i = 0; i < numCells; i++)              for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
475      {                 if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
476        if( elements->Id[i] >= elements->elementDistribution->vtxdist[myRank] && elements->Id[i] <= elements->elementDistribution->vtxdist[myRank+1]-1)                   sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,
477        {                                                    mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
478          sprintf(tmpbuf2,"%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d\n",                                                    mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
479                  nodesGlobal[elements->Nodes[INDEX2(0, i, NN)]],                                                    mesh_p->Nodes->Coordinates[INDEX2(2, i, nDim)]);
480                  nodesGlobal[elements->Nodes[INDEX2(1, i, NN)]],                   __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use);
481                  nodesGlobal[elements->Nodes[INDEX2(2, i, NN)]],                 }
482                  nodesGlobal[elements->Nodes[INDEX2(3, i, NN)]],              }    
483                  nodesGlobal[elements->Nodes[INDEX2(4, i, NN)]],    
484                  nodesGlobal[elements->Nodes[INDEX2(5, i, NN)]],           }
485                  nodesGlobal[elements->Nodes[INDEX2(6, i, NN)]],           #ifdef PASO_MPI
486                  nodesGlobal[elements->Nodes[INDEX2(7, i, NN)]],              if (txt_buffer_in_use==0) { strcpy(txt_buffer, " "); txt_buffer_in_use = 1; } /* avoid zero-length writes */
487                  nodesGlobal[elements->Nodes[INDEX2(8, i, NN)]],              MPI_File_write_ordered(mpi_fileHandle_p, txt_buffer,txt_buffer_in_use, MPI_CHAR, &mpi_status);
488                  nodesGlobal[elements->Nodes[INDEX2(9, i, NN)]],           #endif    
489                  nodesGlobal[elements->Nodes[INDEX2(10, i, NN)]],        } else {
490                  nodesGlobal[elements->Nodes[INDEX2(11, i, NN)]],           if (nDim==2) {
491                  nodesGlobal[elements->Nodes[INDEX2(16, i, NN)]],              for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
492                  nodesGlobal[elements->Nodes[INDEX2(17, i, NN)]],                 if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
493                  nodesGlobal[elements->Nodes[INDEX2(18, i, NN)]],                   fprintf(fileHandle_p,FLOAT_VECTOR_FORMAT,
494                  nodesGlobal[elements->Nodes[INDEX2(19, i, NN)]],                                        mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
495                  nodesGlobal[elements->Nodes[INDEX2(12, i, NN)]],                                        mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
496                  nodesGlobal[elements->Nodes[INDEX2(13, i, NN)]],                                        0.);
497                  nodesGlobal[elements->Nodes[INDEX2(14, i, NN)]],                 }
498                  nodesGlobal[elements->Nodes[INDEX2(15, i, NN)]]);              }      
499          __STRCAT(largebuf,tmpbuf2,tsz)           } else {
500          elementCache.values[pos++]=i;              for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
501        }                 if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
502      }                   fprintf(fileHandle_p,FLOAT_VECTOR_FORMAT,
503    }                                                mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
504    else if (numVTKNodesPerElement!=NN)                                                mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
505    {                                                mesh_p->Nodes->Coordinates[INDEX2(2, i, nDim)]);
506                   }
507      for (i = 0; i < numCells; i++)              }    
508      {    
509        if (elements->Id[i] >= elements->elementDistribution->vtxdist[i] &&  elements->Id[i] <= elements->elementDistribution->vtxdist[i+1] - 1 )           }
510        {        }
511          for (j = 0; j < numVTKNodesPerElement; j++)  
512          {        /* close the Points and open connectivity */
513            sprintf(tmpbuf,"%d ", nodesGlobal[elements->Nodes[INDEX2(elements->ReferenceElement->Type->geoNodes[j], i, NN)]]);  
514            __STRCAT(largebuf,tmpbuf,tsz)        if (mpi_size > 1) {
515          }            if ( my_mpi_rank == 0) {
516          __STRCAT(largebuf,newline,tsz)               #ifdef PASO_MPI
517          elementCache.values[pos++]=i;                  MPI_File_iwrite_shared(mpi_fileHandle_p, tags_End_Points_and_Start_Conn, strlen(tags_End_Points_and_Start_Conn), MPI_CHAR, &mpi_req);
518        }                  MPI_Wait(&mpi_req,&mpi_status);
519      }               #endif
520    }            }
521    else        } else {
522    {           fprintf(fileHandle_p, "%s", tags_End_Points_and_Start_Conn);
523      for(i = 0;i  < numCells ; i++)        }
524      {  
525        if( elements->Id[i] >= elements->elementDistribution->vtxdist[myRank] && elements->Id[i] <= elements->elementDistribution->vtxdist[myRank+1]-1)       /* write the cells */
526        {       if (nodetype == FINLEY_REDUCED_NODES) {
527          for (j = 0; j < numVTKNodesPerElement; j++)          node_index=elements->ReferenceElement->Type->linearNodes;
528          {       } else if (VTK_QUADRATIC_HEXAHEDRON==cellType) {
529            sprintf(tmpbuf,"%d ", nodesGlobal[ elements->Nodes[INDEX2(j, i, NN) ] ] );          node_index=VTK_QUADRATIC_HEXAHEDRON_INDEX;
530            __STRCAT(largebuf,tmpbuf,tsz)       } else if ( (numVTKNodesPerElement!=NN) && (TypeId!=Rec9) && (TypeId!=Hex27) ) {
531            node_index=elements->ReferenceElement->Type->geoNodes;
532         } else {
533            node_index=NULL;
534         }
535    
536         if ( mpi_size > 1) {
537            txt_buffer[0] = '\0';
538            txt_buffer_in_use=0;
539            if (node_index == NULL) {
540               if (TypeId==Rec9) {
541                  for (i = 0; i < numCells; i++) {
542                     if (elements->Owner[i] == my_mpi_rank) {
543                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(0, i, NN)]]);
544                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(4, i, NN)]]);
545                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(8, i, NN)]]);
546                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(7, i, NN)]]);
547                            __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use)
548                            __STRCAT(txt_buffer,NEWLINE,txt_buffer_in_use)
549    
550                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(4, i, NN)]]);
551                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(1, i, NN)]]);
552                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(5, i, NN)]]);
553                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(8, i, NN)]]);
554                            __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use)
555                            __STRCAT(txt_buffer,NEWLINE,txt_buffer_in_use)
556    
557                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(7, i, NN)]]);
558                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(8, i, NN)]]);
559                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(6, i, NN)]]);
560                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(3, i, NN)]]);
561                            __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use)
562                            __STRCAT(txt_buffer,NEWLINE,txt_buffer_in_use)
563    
564                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(8, i, NN)]]);
565                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(5, i, NN)]]);
566                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(2, i, NN)]]);
567                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(6, i, NN)]]);
568                            __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use)
569                            __STRCAT(txt_buffer,NEWLINE,txt_buffer_in_use)
570                    }
571                  }
572               } else if (TypeId==Hex27) {
573                  for (i = 0; i < numCells; i++) {
574                     if (elements->Owner[i] == my_mpi_rank) {
575                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2( 0, i, NN)]]);
576                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2( 8, i, NN)]]);
577                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(20, i, NN)]]);
578                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(11, i, NN)]]);
579                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(12, i, NN)]]);
580                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(21, i, NN)]]);
581                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(26, i, NN)]]);
582                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(24, i, NN)]]);
583                            __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use)
584                            __STRCAT(txt_buffer,NEWLINE,txt_buffer_in_use)
585    
586                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2( 8, i, NN)]]);
587                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2( 1, i, NN)]]);
588                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2( 9, i, NN)]]);
589                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(20, i, NN)]]);
590                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(21, i, NN)]]);
591                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(13, i, NN)]]);
592                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(22, i, NN)]]);
593                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(26, i, NN)]]);
594                            __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use)
595                            __STRCAT(txt_buffer,NEWLINE,txt_buffer_in_use)
596    
597                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(11, i, NN)]]);
598                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(20, i, NN)]]);
599                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(10, i, NN)]]);
600                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2( 3, i, NN)]]);
601                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(24, i, NN)]]);
602                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(26, i, NN)]]);
603                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(23, i, NN)]]);
604                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(15, i, NN)]]);
605                            __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use)
606                            __STRCAT(txt_buffer,NEWLINE,txt_buffer_in_use)
607    
608                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(20, i, NN)]]);
609                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2( 9, i, NN)]]);
610                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2( 2, i, NN)]]);
611                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(10, i, NN)]]);
612                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(26, i, NN)]]);
613                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(22, i, NN)]]);
614                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(14, i, NN)]]);
615                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(23, i, NN)]]);
616                            __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use)
617                            __STRCAT(txt_buffer,NEWLINE,txt_buffer_in_use)
618    
619                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(12, i, NN)]]);
620                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(21, i, NN)]]);
621                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(26, i, NN)]]);
622                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(24, i, NN)]]);
623                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2( 4, i, NN)]]);
624                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(16, i, NN)]]);
625                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(25, i, NN)]]);
626                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(19, i, NN)]]);
627                            __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use)
628                            __STRCAT(txt_buffer,NEWLINE,txt_buffer_in_use)
629    
630                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(21, i, NN)]]);
631                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(13, i, NN)]]);
632                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(22, i, NN)]]);
633                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(26, i, NN)]]);
634                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(16, i, NN)]]);
635                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2( 5, i, NN)]]);
636                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(17, i, NN)]]);
637                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(25, i, NN)]]);
638                            __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use)
639                            __STRCAT(txt_buffer,NEWLINE,txt_buffer_in_use)
640    
641                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(24, i, NN)]]);
642                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(26, i, NN)]]);
643                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(23, i, NN)]]);
644                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(15, i, NN)]]);
645                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(19, i, NN)]]);
646                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(25, i, NN)]]);
647                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(18, i, NN)]]);
648                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2( 7, i, NN)]]);
649                            __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use)
650                            __STRCAT(txt_buffer,NEWLINE,txt_buffer_in_use)
651    
652                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(26, i, NN)]]);
653                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(22, i, NN)]]);
654                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(14, i, NN)]]);
655                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(23, i, NN)]]);
656                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(25, i, NN)]]);
657                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(17, i, NN)]]);
658                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2( 6, i, NN)]]);
659                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(18, i, NN)]]);
660                            __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use)
661                            __STRCAT(txt_buffer,NEWLINE,txt_buffer_in_use)
662                     }
663                  }
664               } else {
665                  for (i = 0; i < numCells; i++) {
666                     if (elements->Owner[i] == my_mpi_rank) {
667                        for (j = 0; j < numVTKNodesPerElement; j++) {
668                            sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(j, i, NN)]]);
669                            __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use)
670                        }
671                        __STRCAT(txt_buffer,NEWLINE,txt_buffer_in_use)
672                     }
673                  }
674               }
675            } else {
676               for (i = 0; i < numCells; i++) {
677                  if (elements->Owner[i] == my_mpi_rank) {
678                     for (j = 0; j < numVTKNodesPerElement; j++) {
679                         sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(node_index[j], i, NN)]]);
680                         __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use)
681                     }
682                     __STRCAT(txt_buffer,NEWLINE,txt_buffer_in_use)
683                  }
684               }
685          }          }
686          __STRCAT(largebuf,newline,tsz)          #ifdef PASO_MPI
687          elementCache.values[pos++]=i;             if (txt_buffer_in_use==0) { strcpy(txt_buffer, " "); txt_buffer_in_use = 1; } /* avoid zero-length writes */
688        }             MPI_File_write_ordered(mpi_fileHandle_p,txt_buffer,txt_buffer_in_use, MPI_CHAR, &mpi_status);
689      }          #endif    
690    }       } else {
691            if (node_index == NULL) {
692    elementCache.size = pos;             if (TypeId==Rec9) {
693                  for (i = 0; i < numCells; i++) {
694    largebuf[tsz] = '\0';                          fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(0, i, NN)]]);
695    MPI_File_write_ordered(fh,largebuf,tsz, MPI_CHAR, &status);                          fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(4, i, NN)]]);
696    MEMFREE(largebuf);                          fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(8, i, NN)]]);
697    MPIO_DEBUG(" Done Writing Connectivity ")                          fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(7, i, NN)]]);
698    MPIO_DEBUG(" Writing Offsets & Types... ")                          fprintf(fileHandle_p,NEWLINE);
699    
700    // Non-Collective                          fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(4, i, NN)]]);
701    if( myRank == 0)                          fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(1, i, NN)]]);
702    {                          fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(5, i, NN)]]);
703      // write out the DataArray element for the offsets                          fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(8, i, NN)]]);
704      char* tag1 = "</DataArray>\n<DataArray Name=\"offsets\" type=\"Int32\" format=\"ascii\">\n";                          fprintf(fileHandle_p,NEWLINE);
705      char* tag2 = "</DataArray>\n";  
706      char *tag3 =  "<DataArray Name=\"types\" type=\"UInt8\" format=\"ascii\">\n";                          fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(7, i, NN)]]);
707      char *tag4 = "</DataArray>\n</Cells>\n";                          fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(8, i, NN)]]);
708                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(6, i, NN)]]);
709      int n = numVTKNodesPerElement;                          fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(3, i, NN)]]);
710                            fprintf(fileHandle_p,NEWLINE);
711      // allocate an upper bound on number of bytes needed    
712      int sz=0;                          fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(8, i, NN)]]);
713      int lg = log10(numGlobalCells * n) + 1;                          fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(5, i, NN)]]);
714      sz += numGlobalCells*lg;                          fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(2, i, NN)]]);
715      sz += numGlobalCells;                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(6, i, NN)]]);
716      tsz = 0;                          fprintf(fileHandle_p,NEWLINE);
717                  }
     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);  
   }  
   
   MPIO_DEBUG(" Done Writing Offsets & Types ")  
718    
719    // Write Cell data header Tags             } else if (TypeId==Hex27) {
720    if(myRank == 0)                   for (i = 0; i < numCells; i++) {
721    {                          fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2( 0, i, NN)]]);
722      MPIO_DEBUG(" Writing Cell Data ...")                          fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2( 8, i, NN)]]);
723      if( write_celldata)                          fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(20, i, NN)]]);
724      {                          fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(11, i, NN)]]);
725        char tmpBuf[80];                          fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(12, i, NN)]]);
726        char header[600];                          fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(21, i, NN)]]);
727        // mark the active data arrays                          fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(26, i, NN)]]);
728        bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;                          fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(24, i, NN)]]);
729        sprintf(tmpBuf, "<CellData");                          fprintf(fileHandle_p,NEWLINE);
730        strcat(header,tmpBuf);  
731        for (i_data =0 ;i_data<num_data;++i_data)                          fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2( 8, i, NN)]]);
732        {                          fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2( 1, i, NN)]]);
733          if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data])                          fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2( 9, i, NN)]]);
734          {                          fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(20, i, NN)]]);
735            // if the rank == 0:   --> scalar data                          fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(21, i, NN)]]);
736            // if the rank == 1:   --> vector data                          fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(13, i, NN)]]);
737            // if the rank == 2:   --> tensor data                          fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(22, i, NN)]]);
738                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(26, i, NN)]]);
739                            fprintf(fileHandle_p,NEWLINE);
740    
741                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(11, i, NN)]]);
742                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(20, i, NN)]]);
743                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(10, i, NN)]]);
744                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2( 3, i, NN)]]);
745                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(24, i, NN)]]);
746                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(26, i, NN)]]);
747                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(23, i, NN)]]);
748                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(15, i, NN)]]);
749                            fprintf(fileHandle_p,NEWLINE);
750    
751                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(20, i, NN)]]);
752                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2( 9, i, NN)]]);
753                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2( 2, i, NN)]]);
754                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(10, i, NN)]]);
755                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(26, i, NN)]]);
756                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(22, i, NN)]]);
757                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(14, i, NN)]]);
758                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(23, i, NN)]]);
759                            fprintf(fileHandle_p,NEWLINE);
760    
761                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(12, i, NN)]]);
762                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(21, i, NN)]]);
763                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(26, i, NN)]]);
764                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(24, i, NN)]]);
765                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2( 4, i, NN)]]);
766                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(16, i, NN)]]);
767                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(25, i, NN)]]);
768                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(19, i, NN)]]);
769                            fprintf(fileHandle_p,NEWLINE);
770    
771                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(21, i, NN)]]);
772                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(13, i, NN)]]);
773                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(22, i, NN)]]);
774                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(26, i, NN)]]);
775                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(16, i, NN)]]);
776                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2( 5, i, NN)]]);
777                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(17, i, NN)]]);
778                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(25, i, NN)]]);
779                            fprintf(fileHandle_p,NEWLINE);
780    
781                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(24, i, NN)]]);
782                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(26, i, NN)]]);
783                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(23, i, NN)]]);
784                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(15, i, NN)]]);
785                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(19, i, NN)]]);
786                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(25, i, NN)]]);
787                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(18, i, NN)]]);
788                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2( 7, i, NN)]]);
789                            fprintf(fileHandle_p,NEWLINE);
790    
791                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(26, i, NN)]]);
792                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(22, i, NN)]]);
793                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(14, i, NN)]]);
794                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(23, i, NN)]]);
795                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(25, i, NN)]]);
796                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(17, i, NN)]]);
797                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2( 6, i, NN)]]);
798                            fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(18, i, NN)]]);
799                            fprintf(fileHandle_p,NEWLINE);
800                     }
801                  } else {
802               for (i = 0; i < numCells; i++) {
803                     for (j = 0; j < numVTKNodesPerElement; j++) {
804                        fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(j, i, NN)]]);
805                      }
806                     fprintf(fileHandle_p,NEWLINE);
807                  }
808              }
809            } else {
810               for (i = 0; i < numCells; i++) {
811                  for (j = 0; j < numVTKNodesPerElement; j++) {
812                     fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(node_index[j], i, NN)]]);
813                   }
814                  fprintf(fileHandle_p,NEWLINE);
815               }
816            }
817    
818         }
819         /* finalize the connection and start the offset section */
820         if (mpi_size > 1) {
821            if( my_mpi_rank == 0) {
822               #ifdef PASO_MPI
823                  MPI_File_iwrite_shared(mpi_fileHandle_p,tags_End_Conn_and_Start_Offset,strlen(tags_End_Conn_and_Start_Offset),MPI_CHAR,&mpi_req);
824                  MPI_Wait(&mpi_req,&mpi_status);
825               #endif
826            }
827         } else {
828            fprintf(fileHandle_p, "%s", tags_End_Conn_and_Start_Offset);
829         }
830    
831        /* write the offsets */
832          
833         if ( mpi_size > 1) {
834            txt_buffer[0] = '\0';
835            txt_buffer_in_use=0;
836            for (i=numVTKNodesPerElement*(myFirstCell*numCellFactor+1); i<=(myFirstCell+myNumCells)*numVTKNodesPerElement*numCellFactor; i+=numVTKNodesPerElement) {
837                sprintf(tmp_buffer, INT_NEWLINE_FORMAT, i);
838                __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use);
839             }
840             #ifdef PASO_MPI
841                if (txt_buffer_in_use==0) { strcpy(txt_buffer, " "); txt_buffer_in_use = 1; } /* avoid zero-length writes */
842                MPI_File_write_ordered(mpi_fileHandle_p,txt_buffer,txt_buffer_in_use, MPI_CHAR, &mpi_status);
843             #endif    
844         } else {
845            for (i=numVTKNodesPerElement; i<=numCells*numVTKNodesPerElement*numCellFactor; i+=numVTKNodesPerElement) {
846               fprintf(fileHandle_p, INT_NEWLINE_FORMAT, i);
847            }
848        
849         }
850         /* finalize the offset section and start the type section */
851         if ( mpi_size > 1) {
852            if ( my_mpi_rank == 0) {
853               #ifdef PASO_MPI
854                  MPI_File_iwrite_shared(mpi_fileHandle_p,tags_End_Offset_and_Start_Type,strlen(tags_End_Offset_and_Start_Type),MPI_CHAR,&mpi_req);
855                  MPI_Wait(&mpi_req,&mpi_status);
856               #endif
857            }
858        } else {
859           fprintf(fileHandle_p, "%s", tags_End_Offset_and_Start_Type);
860        }
861         /* write element type */
862         sprintf(tmp_buffer, INT_NEWLINE_FORMAT, cellType);
863         if ( mpi_size > 1) {
864            txt_buffer[0] = '\0';
865            txt_buffer_in_use=0;
866            for (i=numVTKNodesPerElement*(myFirstCell*numCellFactor+1); i<=(myFirstCell+myNumCells)*numVTKNodesPerElement*numCellFactor; i+=numVTKNodesPerElement) {
867              __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use);
868        }
869             #ifdef PASO_MPI
870                if (txt_buffer_in_use==0) { strcpy(txt_buffer, " "); txt_buffer_in_use = 1; } /* avoid zero-length writes */
871                MPI_File_write_ordered(mpi_fileHandle_p,txt_buffer,txt_buffer_in_use, MPI_CHAR, &mpi_status);
872             #endif    
873         } else {
874            for (i=0; i<numCells*numCellFactor; i++) fprintf(fileHandle_p, "%s",  tmp_buffer);
875         }
876         /* finalize cell information */
877         if ( mpi_size > 1) {
878            if ( my_mpi_rank == 0) {
879               #ifdef PASO_MPI
880                  MPI_File_iwrite_shared(mpi_fileHandle_p,tags_End_Type_And_Cells,strlen(tags_End_Type_And_Cells),MPI_CHAR,&mpi_req);
881                  MPI_Wait(&mpi_req,&mpi_status);
882               #endif
883            }
884        } else {
885           fprintf(fileHandle_p, "%s", tags_End_Type_And_Cells);
886        }
887     }
888    
889     /* Write cell data */
890     if (write_celldata && Finley_noError()) {
891          /* mark the active data arrays */
892          txt_buffer[0] = '\0';
893          set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
894          strcat(txt_buffer, "<CellData");
895          for (i_data =0 ;i_data<num_data;++i_data) {
896            if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data]) {
897              /* if the rank == 0:   --> scalar data */
898              /* if the rank == 1:   --> vector data */
899              /* if the rank == 2:   --> tensor data */
900    
901            switch(getDataPointRank(data_pp[i_data]))            switch(getDataPointRank(data_pp[i_data])) {
           {  
902            case 0:            case 0:
903              if (! set_scalar)              if (! set_scalar) {
904              {                strcat(txt_buffer," Scalars=\"");
905                sprintf(tmpBuf," Scalars=\"%s\"",names_p[i_data]);                strcat(txt_buffer,names_p[i_data]);
906                strcat(header,tmpBuf);                strcat(txt_buffer,"\"");
907                set_scalar=TRUE;                set_scalar=TRUE;
908              }              }
909              break;              break;
910            case 1:            case 1:
911              if (! set_vector)              if (! set_vector) {
912              {                strcat(txt_buffer," Vectors=\"");
913                sprintf(tmpBuf," Vectors=\"%s\"",names_p[i_data]);                strcat(txt_buffer,names_p[i_data]);
914            strcat(header,tmpBuf);                strcat(txt_buffer,"\"");
915                set_vector=TRUE;                set_vector=TRUE;
916              }              }
917              break;              break;
918            case 2:            case 2:
919              if (! set_tensor)              if (! set_tensor) {
920              {                strcat(txt_buffer," Tensors=\"");
921                sprintf(tmpBuf," Tensors=\"%s\"",names_p[i_data]);                strcat(txt_buffer,names_p[i_data]);
922            strcat(header,tmpBuf);                strcat(txt_buffer,"\"");
923                set_tensor=TRUE;                set_tensor=TRUE;
924              }              }
925              break;              break;
# Line 799  void Finley_Mesh_saveVTK_MPIO(const char Line 930  void Finley_Mesh_saveVTK_MPIO(const char
930            }            }
931          }          }
932        }        }
933        strcat(header, ">\n");        strcat(txt_buffer, ">\n");
934        MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);        if ( mpi_size > 1) {
935        MPI_Wait(&req,&status);          if ( my_mpi_rank == 0) {
936      }             #ifdef PASO_MPI
937    }                MPI_File_iwrite_shared(mpi_fileHandle_p,txt_buffer,strlen(txt_buffer),MPI_CHAR,&mpi_req);
938                  MPI_Wait(&mpi_req,&mpi_status);
939    // write actual data (collective)             #endif
940    if(write_celldata)          }
941    {        } else {
942      for (i_data =0 ;i_data<num_data;++i_data)            fprintf(fileHandle_p, "%s", txt_buffer);
943      {        }
944        if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data])        /* write the arrays */
945        {        for (i_data =0 ;i_data<num_data;++i_data) {
946          numPointsPerSample = elements->ReferenceElement->numQuadNodes;           if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data]) {
947          rank = getDataPointRank(data_pp[i_data]);              txt_buffer[0] = '\0';
948          nComp = getDataPointSize(data_pp[i_data]);              txt_buffer_in_use=0;
949          nCompReqd=1;   // the number of components required by vtk              numPointsPerSample=getNumDataPointsPerSample(data_pp[i_data]);
950          shape=0;              rank = getDataPointRank(data_pp[i_data]);
951          if (rank == 0)              nComp = getDataPointSize(data_pp[i_data]);
952          {              nCompReqd=1;   /* the number of components mpi_required by vtk */
953            nCompReqd = 1;              shape=0;
954          }              if (rank == 0) {
955          else if (rank == 1)                nCompReqd = 1;
956          {              } else if (rank == 1) {
957            shape=getDataPointShape(data_pp[i_data], 0);                shape=getDataPointShape(data_pp[i_data], 0);
958            if  (shape>3)                if  (shape>3) {
959            {                  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++;  
960                }                }
961                for (n=0; n<3-shape; n++)                nCompReqd = 3;
962                {              } else {
963                  __STRCAT(largebuf,zero,tsz)                shape=getDataPointShape(data_pp[i_data], 0);
964                  if  (shape>3 || shape != getDataPointShape(data_pp[i_data], 1)) {
965                    Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
966                }                }
967                  nCompReqd = 9;
968              }              }
969              for (m=0; m<3-shape; m++)              if (Finley_noError()) {
970                for (n=0; n<3; n++)                 sprintf(txt_buffer,tag_Float_DataArray,names_p[i_data], nCompReqd);
971                {                 if ( mpi_size > 1) {
972                  __STRCAT(largebuf,zero,tsz)                   if ( my_mpi_rank == 0) {
973                }                      #ifdef PASO_MPI
974            }                         MPI_File_iwrite_shared(mpi_fileHandle_p,txt_buffer,strlen(txt_buffer),MPI_CHAR,&mpi_req);
975            __STRCAT(largebuf,newline,tsz)                         MPI_Wait(&mpi_req,&mpi_status);
976                        #endif
977                     }
978                   } else {
979                       fprintf(fileHandle_p, "%s", txt_buffer);
980                   }
981    
982                   for (i=0; i<numCells; i++) {
983                       if (elements->Owner[i] == my_mpi_rank) {
984                          values = getSampleData(data_pp[i_data], i);
985                          for (l=0; l< numCellFactor;++l) {
986                             /* averaging over the number of points in the sample */
987                             if (isExpanded(data_pp[i_data])) {
988                                  for (k=0; k<MIN(nComp,NCOMP_MAX); k++) sampleAvg[k]=0;
989                                  hits=0;
990                                  for (j=0; j<numPointsPerSample; j++) {
991                                     hits_old=hits;
992                                     if (TypeId==Rec9) {
993                                        switch(l) {
994                                          case 0:
995                                            if (INSIDE_2D(QuadNodes[2*j],QuadNodes[2*j+1],0.25,0.25,0.25)) hits++;  
996                                            break;
997                                          case 1:
998                                            if (INSIDE_2D(QuadNodes[2*j],QuadNodes[2*j+1],0.75,0.25,0.25)) hits++;  
999                                            break;
1000                                          case 2:
1001                                            if (INSIDE_2D(QuadNodes[2*j],QuadNodes[2*j+1],0.25,0.75,0.25)) hits++;  
1002                                            break;
1003                                          case 3:
1004                                            if (INSIDE_2D(QuadNodes[2*j],QuadNodes[2*j+1],0.75,0.75,0.25)) hits++;  
1005                                            break;
1006                                          }
1007                                     } else if (TypeId==Hex27) {
1008                                        switch(l) {
1009                                          case 0:
1010                                            if (INSIDE_3D(QuadNodes[3*j],QuadNodes[3*j+1],QuadNodes[3*j+2],0.25,0.25,0.25,0.25)) hits++;  
1011                                            break;
1012                                          case 1:
1013                                            if (INSIDE_3D(QuadNodes[3*j],QuadNodes[3*j+1],QuadNodes[3*j+2],0.75,0.25,0.25,0.25)) hits++;  
1014                                            break;
1015                                          case 2:
1016                                            if (INSIDE_3D(QuadNodes[3*j],QuadNodes[3*j+1],QuadNodes[3*j+2],0.25,0.75,0.25,0.25)) hits++;  
1017                                            break;
1018                                          case 3:
1019                                            if (INSIDE_3D(QuadNodes[3*j],QuadNodes[3*j+1],QuadNodes[3*j+2],0.75,0.75,0.25,0.25)) hits++;  
1020                                            break;
1021                                          case 4:
1022                                            if (INSIDE_3D(QuadNodes[3*j],QuadNodes[3*j+1],QuadNodes[3*j+2],0.25,0.25,0.75,0.25)) hits++;  
1023                                            break;
1024                                          case 5:
1025                                            if (INSIDE_3D(QuadNodes[3*j],QuadNodes[3*j+1],QuadNodes[3*j+2],0.75,0.25,0.75,0.25)) hits++;  
1026                                            break;
1027                                          case 6:
1028                                            if (INSIDE_3D(QuadNodes[3*j],QuadNodes[3*j+1],QuadNodes[3*j+2],0.25,0.75,0.75,0.25)) hits++;  
1029                                            break;
1030                                          case 7:
1031                                            if (INSIDE_3D(QuadNodes[3*j],QuadNodes[3*j+1],QuadNodes[3*j+2],0.75,0.75,0.75,0.25)) hits++;  
1032                                            break;
1033                                        }
1034                                     } else {
1035                                        hits++;
1036                                     }
1037                                     if (hits_old<hits) for (k=0; k<MIN(nComp,NCOMP_MAX); k++) {
1038                                         sampleAvg[k] += values[INDEX2(k,j,nComp)];
1039                                     }
1040                                  }
1041                                  for (k=0; k<MIN(nComp,NCOMP_MAX); k++) sampleAvg[k] /=MAX(hits,1);
1042                             } else {
1043                                  for (k=0; k<MIN(nComp,NCOMP_MAX); k++) sampleAvg[k] = values[k];
1044                             }
1045                             /* if the number of required components is more than the number
1046                             * of actual components, pad with zeros
1047                             */
1048                             /* probably only need to get shape of first element */
1049                             /* write the data different ways for scalar, vector and tensor */
1050                             if (nCompReqd == 1) {
1051                               sprintf(tmp_buffer,FLOAT_SCALAR_FORMAT,sampleAvg[0]);
1052                             } else if (nCompReqd == 3) {
1053                               if (shape==1) {
1054                                sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,sampleAvg[0],0.,0.);
1055                               } else if (shape==2) {
1056                                sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,sampleAvg[0],sampleAvg[1],0.);
1057                               } else if (shape==3) {
1058                                sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,sampleAvg[0],sampleAvg[1],sampleAvg[2]);
1059                               }
1060                             } else if (nCompReqd == 9) {
1061                               if (shape==1) {
1062                                sprintf(tmp_buffer,FLOAT_TENSOR_FORMAT,sampleAvg[0],0.,0.,
1063                                                                       0.,0.,0.,
1064                                                                    0.,0.,0.);
1065                               } else if (shape==2) {
1066                                sprintf(tmp_buffer,FLOAT_TENSOR_FORMAT,sampleAvg[0],sampleAvg[1],0.,
1067                                                                       sampleAvg[2],sampleAvg[3],0.,
1068                                                                       0.,0.,0.);
1069                               } else if (shape==3) {
1070                                sprintf(tmp_buffer,FLOAT_TENSOR_FORMAT,sampleAvg[0],sampleAvg[1],sampleAvg[2],
1071                                                                       sampleAvg[3],sampleAvg[4],sampleAvg[5],
1072                                                                       sampleAvg[6],sampleAvg[7],sampleAvg[8]);
1073                               }
1074                             }
1075                             /* this needs a bit mor work!!! */
1076                                if ( mpi_size > 1) {
1077                                  __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use);
1078                                } else {
1079                                  fprintf(fileHandle_p, "%s", tmp_buffer);
1080                                }
1081                            }
1082                       }
1083                   }
1084                   if ( mpi_size > 1) {
1085                         #ifdef PASO_MPI
1086                            if (txt_buffer_in_use==0) { strcpy(txt_buffer, " "); txt_buffer_in_use = 1; } /* avoid zero-length writes */
1087                            MPI_File_write_ordered(mpi_fileHandle_p,txt_buffer,txt_buffer_in_use, MPI_CHAR, &mpi_status);
1088                         #endif    
1089                         if ( my_mpi_rank == 0) {
1090                            #ifdef PASO_MPI
1091                               MPI_File_iwrite_shared(mpi_fileHandle_p,tag_End_DataArray,strlen(tag_End_DataArray),MPI_CHAR,&mpi_req);
1092                               MPI_Wait(&mpi_req,&mpi_status);
1093                            #endif
1094                         }
1095                   } else {
1096                       fprintf(fileHandle_p, "%s", tag_End_DataArray);
1097                   }
1098                }
1099             }
1100          }
1101          if ( mpi_size > 1) {
1102            if ( my_mpi_rank == 0) {
1103               #ifdef PASO_MPI
1104                  MPI_File_iwrite_shared(mpi_fileHandle_p,tag_End_CellData,strlen(tag_End_CellData),MPI_CHAR,&mpi_req);
1105                  MPI_Wait(&mpi_req,&mpi_status);
1106               #endif
1107          }          }
1108          largebuf[tsz] = '\0';        } else {
1109          MPI_File_write_ordered(fh,largebuf,tsz,MPI_CHAR,&status);            fprintf(fileHandle_p, "%s", tag_End_CellData);
         MEMFREE(largebuf);  
         if( myRank == 0)  
         {  
           char *tag = "</DataArray>\n";  
           MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);  
           MPI_Wait(&req,&status);  
         }  
   
1110        }        }
     }  
     // 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 ")  
1111    }    }
1112      /* point data */
1113      if (write_pointdata && Finley_noError()) {
1114          /* mark the active data arrays */
1115          set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
1116          txt_buffer[0] = '\0';
1117          strcat(txt_buffer, "<PointData");
1118          for (i_data =0 ;i_data<num_data;++i_data) {
1119            if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data]) {
1120              /* if the rank == 0:   --> scalar data */
1121              /* if the rank == 1:   --> vector data */
1122              /* if the rank == 2:   --> tensor data */
1123    
1124              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]))  
           {  
1125            case 0:            case 0:
1126              if (! set_scalar)              if (! set_scalar) {
1127              {                strcat(txt_buffer," Scalars=\"");
1128                sprintf(tmpBuf," Scalars=\"%s\"",names_p[i_data]);                strcat(txt_buffer,names_p[i_data]);
1129                strcat(header,tmpBuf);                strcat(txt_buffer,"\"");
1130                set_scalar=TRUE;                set_scalar=TRUE;
1131              }              }
1132              break;              break;
1133            case 1:            case 1:
1134              if (! set_vector)              if (! set_vector) {
1135              {                strcat(txt_buffer," Vectors=\"");
1136                sprintf(tmpBuf," Vectors=\"%s\"",names_p[i_data]);                strcat(txt_buffer,names_p[i_data]);
1137                strcat(header,tmpBuf);                strcat(txt_buffer,"\"");
1138                set_vector=TRUE;                set_vector=TRUE;
1139              }              }
1140              break;              break;
1141            case 2:            case 2:
1142              if (! set_tensor)              if (! set_tensor) {
1143              {                strcat(txt_buffer," Tensors=\"");
1144                sprintf(tmpBuf," Tensors=\"%s\"",names_p[i_data]);                strcat(txt_buffer,names_p[i_data]);
1145                strcat(header,tmpBuf);                strcat(txt_buffer,"\"");
1146                set_tensor=TRUE;                set_tensor=TRUE;
1147              }              }
1148              break;              break;
# Line 1000  void Finley_Mesh_saveVTK_MPIO(const char Line 1153  void Finley_Mesh_saveVTK_MPIO(const char
1153            }            }
1154          }          }
1155        }        }
1156        strcat(header, ">\n");        strcat(txt_buffer, ">\n");
1157        MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);        if ( mpi_size > 1) {
1158        MPI_Wait(&req,&status);          if ( my_mpi_rank == 0) {
1159      }             #ifdef PASO_MPI
1160    }                MPI_File_iwrite_shared(mpi_fileHandle_p,txt_buffer,strlen(txt_buffer),MPI_CHAR,&mpi_req);
1161                  MPI_Wait(&mpi_req,&mpi_status);
1162    // write actual data             #endif
1163    if(write_pointdata)          }
1164    {        } else {
1165      for (i_data =0 ;i_data<num_data;++i_data)            fprintf(fileHandle_p, "%s", txt_buffer);
1166      {        }
1167        if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])        /* write the arrays */
1168        {        for (i_data =0 ;i_data<num_data;++i_data) {
1169          numPointsPerSample = elements->ReferenceElement->numQuadNodes;           if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data]) {
1170          rank = getDataPointRank(data_pp[i_data]);              txt_buffer[0] = '\0';
1171          nComp = getDataPointSize(data_pp[i_data]);              txt_buffer_in_use=0;
1172          nCompReqd=1;   // the number of components required by vtk              numPointsPerSample=getNumDataPointsPerSample(data_pp[i_data]);
1173          shape=0;              rank = getDataPointRank(data_pp[i_data]);
1174          if (rank == 0)              nComp = getDataPointSize(data_pp[i_data]);
1175          {              if (getFunctionSpaceType(data_pp[i_data]) == FINLEY_REDUCED_NODES) {
1176            nCompReqd = 1;                 nodeMapping=mesh_p->Nodes->reducedNodesMapping;
1177          }              } else {
1178          else if (rank == 1)                 nodeMapping=mesh_p->Nodes->nodesMapping;
         {  
           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)  
               }  
               for (m=0; m<nCompReqd-shape; m++)  
               {  
                 __STRCAT(largebuf,zero,tsz)  
               }  
1179              }              }
1180              else if (nCompReqd == 9)              nCompReqd=1;   /* the number of components mpi_required by vtk */
1181              {              shape=0;
1182                // tensor data, so have a 3x3 matrix to output as a row              if (rank == 0) {
1183                //  of 9 data points                nCompReqd = 1;
1184                count = 0;              } else if (rank == 1) {
1185                for (m=0; m<shape; m++)                shape=getDataPointShape(data_pp[i_data], 0);
1186                {                if  (shape>3) {
1187                  for (n=0; n<shape; n++)                  Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");
                 {  
                   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)  
                 }  
1188                }                }
1189              }                nCompReqd = 3;
             __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)  
 {  
   char error_msg[LenErrorMsg_MAX];  
   /* if there is no mesh we just return */  
   if (mesh_p==NULL) return;  
   
   int i, j, k, numVTKNodesPerElement,i_data,m, count, n, rank,shape, numPoints, cellType, numCells,  
   nDim, numPointsPerSample, nComp, nCompReqd;  
   
   index_t j2;  
   double* values, rtmp;  
   char elemTypeStr[32];  
   
   /* open the file and check handle */  
   
   FILE * 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 */  
   int nodetype=FINLEY_DEGREES_OF_FREEDOM;  
   int elementtype=FINLEY_UNKNOWN;  
   bool_t isCellCentered[num_data],write_celldata=FALSE,write_pointdata=FALSE;  
   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);  
           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);  
           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);  
           return;  
         }  
         isCellCentered[i_data]=FALSE;  
         break;  
       case FINLEY_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;  
         }  
         isCellCentered[i_data]=TRUE;  
         break;  
       case FINLEY_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);  
           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);  
           return;  
         }  
         isCellCentered[i_data]=TRUE;  
         break;  
       case FINLEY_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);  
           return;  
         }  
         isCellCentered[i_data]=TRUE;  
         break;  
       case FINLEY_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);  
           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);  
         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;  
   Finley_ElementFile* elements=NULL;  
   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);  
     return;  
   }  
   /* map finley element type to VTK element type */  
   numCells = elements->numElements;  
   ElementTypeId TypeId;  
   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);  
     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 */  
   
   int 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);  
           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])  
       {  
         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);  
             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);  
             return;  
           }  
           nCompReqd = 9;  
         }  
         fprintf(fileHandle_p, "<DataArray Name=\"%s\" type=\"Float64\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);  
   
         double sampleAvg[nComp];  
         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<nComp; 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;  
1190              } else {              } else {
1191                 sampleAvg[k] = values[k];                shape=getDataPointShape(data_pp[i_data], 0);
1192              }                if  (shape>3 || shape != getDataPointShape(data_pp[i_data], 1)) {
1193                    Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
           }  
           /* 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)  
           {  
             fprintf(fileHandle_p, " %e", sampleAvg[0]);  
           }  
           else if (nCompReqd == 3)  
           {  
             /* write out the data */  
             for (m=0; m<shape; m++) fprintf(fileHandle_p, " %e", sampleAvg[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", sampleAvg[count]);  
                 count++;  
1194                }                }
1195                for (n=0; n<3-shape; n++) fprintf(fileHandle_p, " %e", 0.);                nCompReqd = 9;
1196              }              }
1197              for (m=0; m<3-shape; m++)              if (Finley_noError()) {
1198                for (n=0; n<3; n++) fprintf(fileHandle_p, " %e", 0.);                 sprintf(txt_buffer,tag_Float_DataArray,names_p[i_data], nCompReqd);
1199            }                 if ( mpi_size > 1) {
1200            fprintf(fileHandle_p, "\n");                   if ( my_mpi_rank == 0) {
1201          }                      #ifdef PASO_MPI
1202          fprintf(fileHandle_p, "</DataArray>\n");                         MPI_File_iwrite_shared(mpi_fileHandle_p,txt_buffer,strlen(txt_buffer),MPI_CHAR,&mpi_req);
1203        }                         MPI_Wait(&mpi_req,&mpi_status);
1204      }                      #endif
1205      fprintf(fileHandle_p, "</CellData>\n");                   }
1206                   } else {
1207                       fprintf(fileHandle_p, "%s", txt_buffer);
1208                   }
1209                   for (i=0; i<mesh_p->Nodes->numNodes; i++) {
1210                      k=globalNodeIndex[i];
1211                      if ( (myFirstNode <= k) && (k < myLastNode) ) {
1212                         values = getSampleData(data_pp[i_data], nodeMapping->target[i]);
1213                         /* if the number of mpi_required components is more than the number
1214                         * of actual components, pad with zeros
1215                         */
1216                         /* probably only need to get shape of first element */
1217                         /* write the data different ways for scalar, vector and tensor */
1218                         if (nCompReqd == 1) {
1219                           sprintf(tmp_buffer,FLOAT_SCALAR_FORMAT,values[0]);
1220                         } else if (nCompReqd == 3) {
1221                           if (shape==1) {
1222                            sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,values[0],0.,0.);
1223                           } else if (shape==2) {
1224                            sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,values[0],values[1],0.);
1225                           } else if (shape==3) {
1226                            sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,values[0],values[1],values[2]);
1227                           }
1228                         } else if (nCompReqd == 9) {
1229                           if (shape==1) {
1230                            sprintf(tmp_buffer,FLOAT_TENSOR_FORMAT,values[0],0.,0.,
1231                                                                   0.,0.,0.,
1232                                                                   0.,0.,0.);
1233                           } else if (shape==2) {
1234                            sprintf(tmp_buffer,FLOAT_TENSOR_FORMAT,values[0],values[1],0.,
1235                                                                   values[2],values[3],0.,
1236                                                                   0.,0.,0.);
1237                           } else if (shape==3) {
1238                            sprintf(tmp_buffer,FLOAT_TENSOR_FORMAT,values[0],values[1],values[2],
1239                                                                   values[3],values[4],values[5],
1240                                                                   values[6],values[7],values[8]);
1241                           }
1242                         }
1243                         if ( mpi_size > 1) {
1244                           __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use);
1245                         } else {
1246                           fprintf(fileHandle_p, "%s", tmp_buffer);
1247                         }
1248                      }
1249                   }
1250                   if ( mpi_size > 1) {
1251                       #ifdef PASO_MPI
1252                         if (txt_buffer_in_use==0) { strcpy(txt_buffer, " "); txt_buffer_in_use = 1; } /* avoid zero-length writes */
1253                         MPI_File_write_ordered(mpi_fileHandle_p,txt_buffer,txt_buffer_in_use, MPI_CHAR, &mpi_status);
1254                       #endif    
1255                       if ( my_mpi_rank == 0) {
1256                          #ifdef PASO_MPI
1257                             MPI_File_iwrite_shared(mpi_fileHandle_p,tag_End_DataArray,strlen(tag_End_DataArray),MPI_CHAR,&mpi_req);
1258                             MPI_Wait(&mpi_req,&mpi_status);
1259                          #endif
1260                       }
1261                   } else {
1262                      fprintf(fileHandle_p, "%s", tag_End_DataArray);
1263                   }
1264                }
1265              }
1266            }
1267            if ( mpi_size > 1) {
1268              if ( my_mpi_rank == 0) {
1269                 #ifdef PASO_MPI
1270                    MPI_File_iwrite_shared(mpi_fileHandle_p,tag_End_PointData,strlen(tag_End_PointData),MPI_CHAR,&mpi_req);
1271                    MPI_Wait(&mpi_req,&mpi_status);
1272                 #endif
1273              }
1274            } else {
1275                fprintf(fileHandle_p, "%s", tag_End_PointData);
1276            }
1277      }
1278      if (Finley_noError()) {
1279         if ( mpi_size > 1) {
1280           if ( my_mpi_rank == 0) {
1281              #ifdef PASO_MPI
1282                 MPI_File_iwrite_shared(mpi_fileHandle_p,footer,strlen(footer),MPI_CHAR,&mpi_req);
1283                 MPI_Wait(&mpi_req,&mpi_status);
1284                 #ifdef MPIO_HINTS
1285                   MPI_Info_free(&mpi_info);
1286                   #undef MPIO_HINTS
1287                 #endif
1288              #endif
1289            }
1290         } else {
1291             fprintf(fileHandle_p, "%s", footer);
1292         }
1293    }    }
1294    /* point data */    if ( mpi_size > 1) {
1295    if (write_pointdata)  #ifdef PASO_MPI
1296    {       MPI_File_close(&mpi_fileHandle_p);
1297      /* mark the active data arrays */  #endif
1298      bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;    } else {
1299      fprintf(fileHandle_p, "<PointData");       fclose(fileHandle_p);
     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);  
           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])  
       {  
         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);  
             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);  
             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  
         */  
         bool_t 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");  
1300    }    }
1301    /* finish off the piece */    TMPMEMFREE(isCellCentered);
1302    fprintf(fileHandle_p, "</Piece>\n");    TMPMEMFREE(txt_buffer);
1303    #else
1304    fprintf(fileHandle_p, "</UnstructuredGrid>\n");    /* Don't kill the job if saveVTK() doesn't work */
1305    /* write the xml footer */    fprintf(stderr, "\n\nsaveVTK warning: VTK is not available\n\n\n");
   fprintf(fileHandle_p, "</VTKFile>\n");  
   /* close the file */  
   fclose(fileHandle_p);  
   return;  
 }  
1306  #endif  #endif
1307    }

Legend:
Removed from v.903  
changed lines
  Added in v.2094

  ViewVC Help
Powered by ViewVC 1.1.26