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

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

  ViewVC Help
Powered by ViewVC 1.1.26