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

Legend:
Removed from v.1030  
changed lines
  Added in v.2548

  ViewVC Help
Powered by ViewVC 1.1.26