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

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

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.26