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

Legend:
Removed from v.1062  
changed lines
  Added in v.2141

  ViewVC Help
Powered by ViewVC 1.1.26