/[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 818 by dhawcroft, Sun Aug 27 11:10:34 2006 UTC revision 2814 by caltinay, Tue Dec 8 01:50:29 2009 UTC
# Line 1  Line 1 
 /*  
 ************************************************************  
 *          Copyright 2006 by ACcESS MNRF                   *  
 *                                                          *  
 *              http://www.access.edu.au                    *  
 *       Primary Business: Queensland, Australia            *  
 *  Licensed under the Open Software License version 3.0    *  
 *     http://www.opensource.org/licenses/osl-3.0.php       *  
 *                                                          *  
 ************************************************************  
 */  
   
   
 /**************************************************************/  
   
 /*   writes data and mesh in a vtk file */  
   
 /**************************************************************/  
   
 /*   Author: Paul Cochrane, cochrane@esscc.uq.edu.au */  
 /*   MPI-IO version: Derick Hawcroft, d.hawcroft@uq.edu.au         */  
   
 /*   Version: $Id$ */  
   
 /**************************************************************/  
1    
2    /*******************************************************
3    *
4    * Copyright (c) 2003-2009 by University of Queensland
5    * Earth Systems Science Computational Center (ESSCC)
6    * http://www.uq.edu.au/esscc
7    *
8    * Primary Business: Queensland, Australia
9    * Licensed under the Open Software License version 3.0
10    * http://www.opensource.org/licenses/osl-3.0.php
11    *
12    *******************************************************/
13    
14    /***************************************************************************/
15    /*   Writes data and mesh in VTK XML format to a VTU file.                 */
16    /*   Nodal data needs to be given on FINLEY_NODES or FINLEY_REDUCED_NODES  */
17    /***************************************************************************/
18    
19  #include "Mesh.h"  #include "Mesh.h"
20  #include "vtkCellType.h"  /* copied from vtk source directory !!! */  #include "Assemble.h"
21    #include "vtkCellType.h"  /* copied from vtk source directory */
22  /*  #include "paso/PasoUtil.h"
23   MPI version notes:  
24    #define INT_FORMAT "%d "
25   ******************************************************************************  #define LEN_INT_FORMAT (unsigned int)(9+2)
26   ***                                                                       ****  #define INT_NEWLINE_FORMAT "%d\n"
27   *** WARNING: Won't work for meshes with peridodic boundary conditions yet ****  #define SCALAR_FORMAT "%12.6e\n"
28   ***                                                                       ****    #define VECTOR_FORMAT "%12.6e %12.6e %12.6e\n"
29   ******************************************************************************  #define TENSOR_FORMAT "%12.6e %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e\n"
30    /* strlen("-1.234567e+789 ") == 15 */
31   In this version, the rank==0 process writes *all* opening and closing  #define LEN_TENSOR_FORMAT (unsigned int)(9*15+2)
32   XML tags.  #define NEWLINE "\n"
33   Individual process data is copied to a buffer before being written  #define LEN_TMP_BUFFER LEN_TENSOR_FORMAT+(MAX_numNodes*LEN_INT_FORMAT+1)+2
34   out. The  routines are collectively called and will be called in the natural  #define NCOMP_MAX (unsigned int)9
35   ordering i.e 0 to maxProcs-1.  
36    #define __STRCAT(dest, chunk, dest_in_use) \
37   Notable Notables:  do {\
38   the struct localIndexCache: stores local domain indices for faster  reference      strcpy(&dest[dest_in_use], chunk);\
39  */      dest_in_use += strlen(chunk);\
40    } while(0)
41    
42  #ifdef PASO_MPI  #ifdef PASO_MPI
43    /* writes buffer to file catching the empty buffer case which causes problems
44     * with some MPI versions */
45    #define MPI_WRITE_ORDERED(BUF) \
46    do {\
47        int LLEN=0; \
48        LLEN=(int) strlen(BUF); \
49        if (LLEN==0) { strcpy(BUF, ""); LLEN=0; }\
50        MPI_File_write_ordered(mpi_fileHandle_p, BUF, LLEN, MPI_CHAR, &mpi_status);\
51    } while(0)
52    
53    /* writes buffer to file on master only */
54    #define MPI_RANK0_WRITE_SHARED(BUF) \
55    do {\
56        int LLEN=0; \
57        if (my_mpi_rank == 0) {\
58            LLEN=(int) strlen(BUF); \
59            if (LLEN==0) { strcpy(BUF,""); LLEN=0; }\
60            MPI_File_iwrite_shared(mpi_fileHandle_p, BUF, LLEN, MPI_CHAR, &mpi_req);\
61            MPI_Wait(&mpi_req, &mpi_status);\
62        }\
63    } while(0)
64    
65    /* For reference only. Investigation needed as to which values may improve
66     * performance */
67    #if 0
68    void create_MPIInfo(MPI_Info& info)
69    {
70        MPI_Info_create(&info);
71        MPI_Info_set(info, "access_style", "write_once, sequential");
72        MPI_Info_set(info, "collective_buffering", "true");
73        MPI_Info_set(info, "cb_block_size",        "131072");
74        MPI_Info_set(info, "cb_buffer_size",       "1048567");
75        MPI_Info_set(info, "cb_nodes",             "8");
76        MPI_Info_set(info, "striping_factor",      "16");
77        MPI_Info_set(info, "striping_unit",        "424288");
78    }
79    #endif
80    
81    #else
82    
83  //#define MPIO_HINTS  #define MPI_WRITE_ORDERED(A)
84    #define MPI_RANK0_WRITE_SHARED(A)
85    
86    #endif /* PASO_MPI */
87    
 #define MPIO_DEBUG(str) \  
 { \  
     if(myRank == 0) \  
     printf("==== MPI-IO => %s \n", str); \  
 }  
88    
89  void Finley_Mesh_saveVTK_MPIO(const char * filename_p, Finley_Mesh *mesh_p, const dim_t num_data,char* *names_p, escriptDataC* *data_pp)  /* Returns one if the node given by coords and idx is within the quadrant
90     * indexed by q and if the element type is Rec9 or Hex27, zero otherwise */
91    int nodeInQuadrant(const double *coords, ElementTypeId type, int idx, int q)
92  {  {
93    int    numPoints,  #define INSIDE_1D(_X_,_C_,_R_) ( ABS((_X_)-(_C_)) <= (_R_) )
94    numCells = -1,  #define INSIDE_2D(_X_,_Y_,_CX_,_CY_,_R_) ( INSIDE_1D(_X_,_CX_,_R_) && INSIDE_1D(_Y_,_CY_,_R_))
95               myRank,comm,gsize,  #define INSIDE_3D(_X_,_Y_,_Z_,_CX_,_CY_,_CZ_,_R_) ( INSIDE_1D(_X_,_CX_,_R_) && INSIDE_1D(_Y_,_CY_,_R_) && INSIDE_1D(_Z_,_CZ_,_R_) )
96               numLocal,  
97               nDim,      int ret;
98               shape;      if ( type == 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;  
         }  
         else  
         {  
           Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");  
           return ;  
         }  
         isCellCentered[i_data]=FALSE;  
         break;  
       case FINLEY_REDUCED_DEGREES_OF_FREEDOM:  
         nodetype = FINLEY_REDUCED_DEGREES_OF_FREEDOM;  
         if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)  
         {  
           elementtype=FINLEY_ELEMENTS;  
         }  
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_NODES:          else if (q==2)
111          nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_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_ELEMENTS:              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          else              ret = INSIDE_3D(coords[3*idx], coords[3*idx+1], coords[3*idx+2],
128          {                      0.75, 0.75, 0.25, 0.25);
129            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");          else if (q==4)
130            return;              ret = INSIDE_3D(coords[3*idx], coords[3*idx+1], coords[3*idx+2],
131          }                      0.25, 0.25, 0.75, 0.25);
132          isCellCentered[i_data]=TRUE;          else if (q==5)
133          break;              ret = INSIDE_3D(coords[3*idx], coords[3*idx+1], coords[3*idx+2],
134        case FINLEY_FACE_ELEMENTS:                      0.75, 0.25, 0.75, 0.25);
135          nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;          else if (q==6)
136          if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_FACE_ELEMENTS)              ret = INSIDE_3D(coords[3*idx], coords[3*idx+1], coords[3*idx+2],
137          {                      0.25, 0.75, 0.75, 0.25);
138            elementtype=FINLEY_FACE_ELEMENTS;          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        return ret;
147    }
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=\"Float32\" format=\"ascii\">\n"  
             ,numPoints,numGlobalCells,MAX(3,nDim));  
   
   
     MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);  
     MPI_Wait(&req,&status);  
   }  
   
   MPIO_DEBUG(" Writing Coordinate Points... ")  
   
   numLocalNodes=localDOF;  
   
   //  we will be writing values which vary from 13-15 chars hence the strlen()  
   char* largebuf = MEMALLOC( numLocalNodes*15*nDim + numLocalNodes*2 + 1 ,char);  
   largebuf[0] = '\0';  
   char tmpbuf[15];  
   int tsz=0;  
   int numNodesOutput=0;  
   index_t pos=0;  
   
   index_t *vtxdist = NULL, *DOFNodes=NULL,*forwardBuffer=NULL,*backwardBuffer=NULL;  
   
   DOFNodes   = MEMALLOC(mesh_p->Nodes->numNodes,index_t);  
   nodeCache.values = MEMALLOC( numLocalNodes, index_t);  
   index_t bc_pos = 0;  
   
   // Custom string concat:  avoids expensive strlen(3) call by strcat(3)  
   // Note the implicit assumption on the variable "tsz"  
   int __len,__j;  
   char  *zero = "0.000000e+00 ";  
   char  *newline = "\n";  
     
 #define __STRCAT(dest,chunk,tsz)  \  
 {                  \  
    __len = strlen(chunk); \  
    __j = -1;      \  
    while(__j++ < __len)  \  
     *(dest+tsz+__j)=*(chunk+__j); \  
    tsz+=__len;              \  
 }  
     
   // Loop over all nodes      
   for (i = 0; i < mesh_p->Nodes->numNodes; i++)  
   {  
     // This is the bit that will break for periodic BCs because it assumes that there is a one to one  
     // correspondance between nodes and Degrees of freedom  
     //TODO: handle periodic BC's  
     DOFNodes[mesh_p->Nodes->degreeOfFreedom[i]] = i;  
   
     // Is this node local to the domain ?  
     if( mesh_p->Nodes->degreeOfFreedom[i] < localDOF )  
     {  
       for (j = 0; j < nDim; j++)  
       {  
         sprintf(tmpbuf,"%e ", mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)] );  
         __STRCAT(largebuf,tmpbuf,tsz)  
       }  
       for (k=0; k<3-nDim; k++)  
       {  
         __STRCAT(largebuf,zero,tsz)  
       }  
       __STRCAT(largebuf,newline,tsz)  
       nodeCache.values[numNodesOutput++]=i;  
     }  
   }  
   
   nodeCache.size=numNodesOutput;  
   
   largebuf[tsz] = '\0';  
   MPI_File_write_ordered(fh, largebuf,tsz, MPI_CHAR, &status);  
   MEMFREE(largebuf);  
   
   nodesGlobal = MEMALLOC(mesh_p->Nodes->numNodes ,index_t);  
   
   // form distribution info on who output which nodes  
   vtxdist = MEMALLOC( gsize+1, index_t );  
   vtxdist[0]=0;  
   MPI_Allgather(&numNodesOutput,1,MPI_INT,vtxdist+1,1,MPI_INT,comm);  
   for( i=0; i<gsize; i++ )  
     vtxdist[i+1]+=vtxdist[i];  
   
   // will not work for periodic boundary conditions  
   // calculate the local nodes file positions  
   pos = 0;  
   for( i=0; i<mesh_p->Nodes->numNodes; i++ )  
   {  
     if( mesh_p->Nodes->degreeOfFreedom[i]< localDOF )  
     {  
       nodesGlobal[i] = vtxdist[myRank] + pos++;  
     }  
     else  
       nodesGlobal[i] = -1;  
   }  
   
   // communicate the local Nodes file position to the interested parties  
   // send local info  
   forwardBuffer = MEMALLOC( mesh_p->Nodes->numNodes, index_t );  
   for( n=0; n < dist->numNeighbours; n++ )  
   {  
     if(  dist->edges[n]->numForward)  
     {  
       for( i=0; i < dist->edges[n]->numForward; i++ )  
         forwardBuffer[i] = nodesGlobal[DOFNodes[dist->edges[n]->indexForward[i] ]];  
       Paso_CommBuffer_pack( mesh_p->Nodes->CommBuffer, dist->neighbours[n], NULL, forwardBuffer, sizeof(index_t), 0 );  
       Paso_CommBuffer_send( mesh_p->Nodes->CommBuffer, dist->neighbours[n], sizeof(index_t) );  
     }  
   }  
   // receive external info  
   backwardBuffer = MEMALLOC( mesh_p->Nodes->numNodes, index_t );  
   for( n=0; n < dist->numNeighbours; n++ )  
   {  
     if( dist->edges[n]->numBackward )  
     {  
       Paso_CommBuffer_recv(mesh_p->Nodes->CommBuffer, dist->neighbours[n], sizeof(index_t));  
       Paso_CommBuffer_unpack(mesh_p->Nodes->CommBuffer, dist->neighbours[n], NULL, backwardBuffer, sizeof(index_t), 0 );  
       /* TODO: voodoo to handle perdiodic  BC's */  
       for( i=0; i<dist->edges[n]->numBackward; i++ )  
         nodesGlobal[DOFNodes[dist->edges[n]->indexBackward[i] ]] = backwardBuffer[i];  
221      }      }
222    }      my_mpi_rank = mesh_p->Nodes->MPIInfo->rank;
223          mpi_size = mesh_p->Nodes->MPIInfo->size;
224    
225          /************************************************************************
226    MEMFREE(vtxdist);       * open the file and check handle *
227    MEMFREE(DOFNodes);       */
228    MEMFREE(backwardBuffer);      if (mpi_size > 1) {
229    MEMFREE(forwardBuffer);  #ifdef PASO_MPI
230            const int amode = MPI_MODE_CREATE|MPI_MODE_WRONLY|MPI_MODE_UNIQUE_OPEN;
231    if( myRank == 0)          int ierr;
232    {          if (my_mpi_rank == 0 && Paso_fileExists(filename_p)) {
233      char* tags = "</DataArray>\n</Points>\n<Cells>\n<DataArray Name=\"connectivity\" type=\"Int32\" " \              remove(filename_p);
234                   "format=\"ascii\">\n" ;          }
235      MPI_File_iwrite_shared(fh,tags,strlen(tags),MPI_CHAR,&req);          ierr = MPI_File_open(mesh_p->Nodes->MPIInfo->comm, (char*)filename_p,
236      MPI_Wait(&req,&status);                               amode, mpi_info, &mpi_fileHandle_p);
237    }          if (ierr != MPI_SUCCESS) {
238    MPIO_DEBUG(" Done Writing Coordinate Points ")              sprintf(errorMsg, "saveVTK: File %s could not be opened for writing in parallel.", filename_p);
239                Finley_setError(IO_ERROR, errorMsg);
240    /* BEGIN CONNECTIVITY */          } else {
241                ierr=MPI_File_set_view(mpi_fileHandle_p,MPI_DISPLACEMENT_CURRENT,
242    int NN = elements->ReferenceElement->Type->numNodes; /* num Nodes holding ref-element */                      MPI_CHAR, MPI_CHAR, "native", mpi_info);
243            }
244    // Collective  #endif /* PASO_MPI */
245    MPIO_DEBUG(" Writing Connectivity... ")      } else {
246            fileHandle_p = fopen(filename_p, "w");
247    size_t sz = numLocalCells*6*numVTKNodesPerElement + numLocalCells;          if (fileHandle_p==NULL) {
248    largebuf = MEMALLOC(sz,char);              sprintf(errorMsg, "saveVTK: File %s could not be opened for writing.", filename_p);
249    largebuf[0] = '\0';              Finley_setError(IO_ERROR, errorMsg);
250    tsz=0;          }
251    pos = 0;      }
252    // numCells?      if (!Paso_MPIInfo_noError(mesh_p->Nodes->MPIInfo)) return;
253    elementCache.values = MEMALLOC(numLocalCells,index_t);  
254    if (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM)      /* General note: From this point if an error occurs Finley_setError is
255    {       * called and subsequent steps are skipped until the end of this function
256      for (i = 0; i < numCells; i++)       * where allocated memory is freed and the file is closed. */
257      {  
258        /************************************************************************/
259        if (elements->Id[i] >= elements->elementDistribution->vtxdist[i] &&  elements->Id[i] <= elements->elementDistribution->vtxdist[i+1] - 1 )      /* find the mesh type to be written */
260        {  
261          for (j = 0; j < numVTKNodesPerElement; j++)      isCellCentered = TMPMEMALLOC(num_data, bool_t);
262          {      maxNameLen = 0;
263            sprintf(tmpbuf,"%d ",nodesGlobal[mesh_p->Nodes->toReduced[elements->Nodes[INDEX2(elements->ReferenceElement->Type->linearNodes[j], i, NN)]]]);      if (!Finley_checkPtr(isCellCentered)) {
264            __STRCAT(largebuf,tmpbuf,tsz)          for (dataIdx=0; dataIdx<num_data; ++dataIdx) {
265          }              if (! isEmpty(data_pp[dataIdx])) {
266          __STRCAT(largebuf,newline,tsz)                  switch(getFunctionSpaceType(data_pp[dataIdx]) ) {
267          elementCache.values[pos++]=i;                      case FINLEY_NODES:
268        }                          nodeType = (nodeType == FINLEY_REDUCED_NODES) ? FINLEY_REDUCED_NODES : FINLEY_NODES;
269      }                          isCellCentered[dataIdx] = FALSE;
270    }                          if (elementType==FINLEY_UNKNOWN || elementType==FINLEY_ELEMENTS) {
271    else if (VTK_QUADRATIC_HEXAHEDRON==cellType)                              elementType = FINLEY_ELEMENTS;
272    {                          } else {
273      char tmpbuf2[20*20*2];                              Finley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
274                            }
275      for (i = 0; i < numCells; i++)                      break;
276      {                      case FINLEY_REDUCED_NODES:
277        if( elements->Id[i] >= elements->elementDistribution->vtxdist[myRank] && elements->Id[i] <= elements->elementDistribution->vtxdist[myRank+1]-1)                          nodeType = FINLEY_REDUCED_NODES;
278        {                          isCellCentered[dataIdx] = FALSE;
279          sprintf(tmpbuf2,"%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d\n",                          if (elementType==FINLEY_UNKNOWN || elementType==FINLEY_ELEMENTS) {
280                  nodesGlobal[elements->Nodes[INDEX2(0, i, NN)]],                              elementType = FINLEY_ELEMENTS;
281                  nodesGlobal[elements->Nodes[INDEX2(1, i, NN)]],                          } else {
282                  nodesGlobal[elements->Nodes[INDEX2(2, i, NN)]],                              Finley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
283                  nodesGlobal[elements->Nodes[INDEX2(3, i, NN)]],                      }
284                  nodesGlobal[elements->Nodes[INDEX2(4, i, NN)]],                      break;
285                  nodesGlobal[elements->Nodes[INDEX2(5, i, NN)]],                      case FINLEY_REDUCED_ELEMENTS:
286                  nodesGlobal[elements->Nodes[INDEX2(6, i, NN)]],                          hasReducedElements = TRUE;
287                  nodesGlobal[elements->Nodes[INDEX2(7, i, NN)]],                      case FINLEY_ELEMENTS:
288                  nodesGlobal[elements->Nodes[INDEX2(8, i, NN)]],                          isCellCentered[dataIdx] = TRUE;
289                  nodesGlobal[elements->Nodes[INDEX2(9, i, NN)]],                          if (elementType==FINLEY_UNKNOWN || elementType==FINLEY_ELEMENTS) {
290                  nodesGlobal[elements->Nodes[INDEX2(10, i, NN)]],                              elementType = FINLEY_ELEMENTS;
291                  nodesGlobal[elements->Nodes[INDEX2(11, i, NN)]],                          } else {
292                  nodesGlobal[elements->Nodes[INDEX2(16, i, NN)]],                              Finley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
293                  nodesGlobal[elements->Nodes[INDEX2(17, i, NN)]],                          }
294                  nodesGlobal[elements->Nodes[INDEX2(18, i, NN)]],                      break;
295                  nodesGlobal[elements->Nodes[INDEX2(19, i, NN)]],                      case FINLEY_REDUCED_FACE_ELEMENTS:
296                  nodesGlobal[elements->Nodes[INDEX2(12, i, NN)]],                          hasReducedElements = TRUE;
297                  nodesGlobal[elements->Nodes[INDEX2(13, i, NN)]],                      case FINLEY_FACE_ELEMENTS:
298                  nodesGlobal[elements->Nodes[INDEX2(14, i, NN)]],                          isCellCentered[dataIdx] = TRUE;
299                  nodesGlobal[elements->Nodes[INDEX2(15, i, NN)]]);                          if (elementType==FINLEY_UNKNOWN || elementType==FINLEY_FACE_ELEMENTS) {
300          __STRCAT(largebuf,tmpbuf2,tsz)                              elementType = FINLEY_FACE_ELEMENTS;
301          elementCache.values[pos++]=i;                          } else {
302        }                              Finley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
303      }                          }
304    }                      break;
305    else if (numVTKNodesPerElement!=NN)                      case FINLEY_POINTS:
306    {                          isCellCentered[dataIdx]=TRUE;
307                            if (elementType==FINLEY_UNKNOWN || elementType==FINLEY_POINTS) {
308      for (i = 0; i < numCells; i++)                              elementType = FINLEY_POINTS;
309      {                          } else {
310        if (elements->Id[i] >= elements->elementDistribution->vtxdist[i] &&  elements->Id[i] <= elements->elementDistribution->vtxdist[i+1] - 1 )                              Finley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
311        {                          }
312          for (j = 0; j < numVTKNodesPerElement; j++)                      break;
313          {                      case FINLEY_REDUCED_CONTACT_ELEMENTS_1:
314            sprintf(tmpbuf,"%d ", nodesGlobal[elements->Nodes[INDEX2(elements->ReferenceElement->Type->geoNodes[j], i, NN)]]);                          hasReducedElements = TRUE;
315            __STRCAT(largebuf,tmpbuf,tsz)                      case FINLEY_CONTACT_ELEMENTS_1:
316          }                          isCellCentered[dataIdx] = TRUE;
317          __STRCAT(largebuf,newline,tsz)                          if (elementType==FINLEY_UNKNOWN || elementType==FINLEY_CONTACT_ELEMENTS_1) {
318          elementCache.values[pos++]=i;                              elementType = FINLEY_CONTACT_ELEMENTS_1;
319        }                          } else {
320      }                              Finley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
321    }                          }
322    else                      break;
323    {                      case FINLEY_REDUCED_CONTACT_ELEMENTS_2:
324      for(i = 0;i  < numCells ; i++)                          hasReducedElements = TRUE;
325      {                      case FINLEY_CONTACT_ELEMENTS_2:
326        if( elements->Id[i] >= elements->elementDistribution->vtxdist[myRank] && elements->Id[i] <= elements->elementDistribution->vtxdist[myRank+1]-1)                          isCellCentered[dataIdx] = TRUE;
327        {                          if (elementType==FINLEY_UNKNOWN || elementType==FINLEY_CONTACT_ELEMENTS_1) {
328          for (j = 0; j < numVTKNodesPerElement; j++)                              elementType = FINLEY_CONTACT_ELEMENTS_1;
329          {                          } else {
330            sprintf(tmpbuf,"%d ", nodesGlobal[ elements->Nodes[INDEX2(j, i, NN) ] ] );                              Finley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
331            __STRCAT(largebuf,tmpbuf,tsz)                          }
332                        break;
333                        default:
334                            sprintf(errorMsg, "saveVTK: unknown function space type %d",getFunctionSpaceType(data_pp[dataIdx]));
335                            Finley_setError(TYPE_ERROR, errorMsg);
336                    }
337                    if (isCellCentered[dataIdx]) {
338                        writeCellData = TRUE;
339                    } else {
340                        writePointData = TRUE;
341                    }
342                    maxNameLen = MAX(maxNameLen, strlen(names_p[dataIdx]));
343                }
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;  
           }  
377          }          }
378        }          if (elements==NULL) {
379        strcat(header, ">\n");              Finley_setError(SYSTEM_ERROR, "saveVTK: undefined element file");
380        MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);          } else {
381        MPI_Wait(&req,&status);              /* map finley element type to VTK element type */
382      }              numCells = elements->numElements;
383    }              globalNumCells = Finley_ElementFile_getGlobalNumElements(elements);
384                myNumCells = Finley_ElementFile_getMyNumElements(elements);
385                myFirstCell = Finley_ElementFile_getFirstElement(elements);
386                NN = elements->numNodes;
387                if (hasReducedElements) {
388                    quadNodes_p=elements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->QuadNodes;
389                } else {
390                    quadNodes_p=elements->referenceElementSet->referenceElement->Parametrization->QuadNodes;
391                }
392                if (nodeType==FINLEY_REDUCED_NODES) {
393                    typeId = elements->referenceElementSet->referenceElement->LinearType->TypeId;
394                } else {
395                    typeId = elements->referenceElementSet->referenceElement->Type->TypeId;
396                }
397                switch (typeId) {
398                    case Point1:
399                    case Line2Face:
400                    case Line3Face:
401                    case Point1_Contact:
402                    case Line2Face_Contact:
403                    case Line3Face_Contact:
404                        cellType = VTK_VERTEX;
405                        numVTKNodesPerElement = 1;
406                    break;
407    
408    // write actual data (collective)                  case Line2:
409    if(write_celldata)                  case Tri3Face:
410    {                  case Rec4Face:
411      for (i_data =0 ;i_data<num_data;++i_data)                  case Line2_Contact:
412      {                  case Tri3_Contact:
413        if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data])                  case Tri3Face_Contact:
414        {                  case Rec4Face_Contact:
415          numPointsPerSample = elements->ReferenceElement->numQuadNodes;                      cellType = VTK_LINE;
416          rank = getDataPointRank(data_pp[i_data]);                      numVTKNodesPerElement = 2;
417          nComp = getDataPointSize(data_pp[i_data]);                  break;
         nCompReqd=1;   // the number of components required by vtk  
         shape=0;  
         if (rank == 0)  
         {  
           nCompReqd = 1;  
         }  
         else if (rank == 1)  
         {  
           shape=getDataPointShape(data_pp[i_data], 0);  
           if  (shape>3)  
           {  
             Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");  
             return;  
           }  
           nCompReqd = 3;  
         }  
         else  
         {  
           shape=getDataPointShape(data_pp[i_data], 0);  
           if  (shape>3 || shape != getDataPointShape(data_pp[i_data], 1))  
           {  
             Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");  
             return;  
           }  
           nCompReqd = 9;  
         }  
418    
419          if( myRank == 0)          case Line3Macro:
420          {                      cellType = VTK_LINE;
421            char header[250];                      numCellFactor = 2;
422            sprintf(header,"<DataArray Name=\"%s\" type=\"Float32\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);                      numVTKNodesPerElement = 2;
423            MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);                  break;
           MPI_Wait(&req,&status);  
         }  
424    
425          // Write the actual data                  case Tri3:
426          char tmpbuf[15];                  case Tet4Face:
427          char* largebuf = MEMALLOC(nCompReqd*numLocalCells*15 + numLocalCells*nCompReqd + nCompReqd + 15,char);                  case Tet4Face_Contact:
428          largebuf[0] = '\0';                      cellType = VTK_TRIANGLE;
429          size_t tsz = 0;                      numVTKNodesPerElement = 3;
430                    break;
         double sampleAvg[nComp];  
   
         for (k=0; k<elementCache.size; k++)  
         {  
           i = elementCache.values[k];  
   
           values = getSampleData(data_pp[i_data], i);  
           // averaging over the number of points in the sample  
           for (n=0; n<nComp; n++)  
           {  
             rtmp = 0.;  
             for (j=0; j<numPointsPerSample; j++) rtmp += values[INDEX2(n,j,nComp)];  
             sampleAvg[k] = rtmp/numPointsPerSample;  
           }  
           // if the number of required components is more than the number  
           // of actual components, pad with zeros  
   
           // probably only need to get shape of first element  
           // write the data different ways for scalar, vector and tensor  
           if (nCompReqd == 1)  
           {  
             sprintf(tmpbuf, " %e", sampleAvg[0]);  
             __STRCAT(largebuf,tmpbuf,tsz)  
           }  
           else if (nCompReqd == 3)  
           {  
             // write out the data  
             for (m=0; m<shape; m++)  
             {  
               sprintf(tmpbuf, " %e", sampleAvg[m]);  
               __STRCAT(largebuf,tmpbuf,tsz)  
             }  
             for (m=0; m<nCompReqd-shape; m++)  
             {  
               __STRCAT(largebuf,zero,tsz)  
             }  
           }  
           else if (nCompReqd == 9)  
           {  
             // tensor data, so have a 3x3 matrix to output as a row  
             // of 9 data points  
             count = 0;  
             for (m=0; m<shape; m++)  
             {  
               for (n=0; n<shape; n++)  
               {  
                 sprintf(tmpbuf, " %e", sampleAvg[count]);  
                 __STRCAT(largebuf,tmpbuf,tsz)  
                 count++;  
               }  
               for (n=0; n<3-shape; n++)  
               {  
                 __STRCAT(largebuf,zero,tsz)  
               }  
             }  
             for (m=0; m<3-shape; m++)  
               for (n=0; n<3; n++)  
               {  
                 __STRCAT(largebuf,zero,tsz)  
               }  
           }  
           __STRCAT(largebuf,newline,tsz)  
         }  
         largebuf[tsz] = '\0';  
         MPI_File_write_ordered(fh,largebuf,tsz,MPI_CHAR,&status);  
         MEMFREE(largebuf);  
         if( myRank == 0)  
         {  
           char *tag = "</DataArray>\n";  
           MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);  
           MPI_Wait(&req,&status);  
         }  
431    
432        }                  case Rec4:
433      }                  case Hex8Face:
434      // closing celldata tag                  case Rec4_Contact:
435      if(myRank == 0)                  case Hex8Face_Contact:
436      {                      cellType = VTK_QUAD;
437        char* tag =  "</CellData>\n";                      numVTKNodesPerElement = 4;
438        MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);                  break;
       MPI_Wait(&req,&status);  
     }  
439    
440      MPIO_DEBUG(" Done Writing Cell Data ")                  case Rec9Macro:
441    }          case Rec9:
442                        numCellFactor = 4;
443                        cellType = VTK_QUAD;
444                        numVTKNodesPerElement = 4;
445                    break;
446    
447                    case Tet4:
448                        cellType = VTK_TETRA;
449                        numVTKNodesPerElement = 4;
450                    break;
451    
452    // Write Point Data Header Tags                  case Hex8:
453    if( myRank == 0)                      cellType = VTK_HEXAHEDRON;
454    {                      numVTKNodesPerElement = 8;
455      char header[600];                  break;
     char tmpBuf[50];  
   
     if (write_pointdata)  
     {  
       MPIO_DEBUG(" Writing Pointdata... ")  
       // mark the active data arrays  
       bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;  
       sprintf(header, "<PointData");  
       for (i_data =0 ;i_data<num_data;++i_data)  
       {  
         if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])  
         {  
           // if the rank == 0:   --> scalar data  
           // if the rank == 1:   --> vector data  
           // if the rank == 2:   --> tensor data  
   
           switch(getDataPointRank(data_pp[i_data]))  
           {  
           case 0:  
             if (! set_scalar)  
             {  
               sprintf(tmpBuf," Scalars=\"%s\"",names_p[i_data]);  
               strcat(header,tmpBuf);  
               set_scalar=TRUE;  
             }  
             break;  
           case 1:  
             if (! set_vector)  
             {  
               sprintf(tmpBuf," Vectors=\"%s\"",names_p[i_data]);  
               strcat(header,tmpBuf);  
               set_vector=TRUE;  
             }  
             break;  
           case 2:  
             if (! set_tensor)  
             {  
               sprintf(tmpBuf," Tensors=\"%s\"",names_p[i_data]);  
               strcat(header,tmpBuf);  
               set_tensor=TRUE;  
             }  
             break;  
           default:  
             sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);  
             Finley_setError(VALUE_ERROR,error_msg);  
             return;  
           }  
         }  
       }  
       strcat(header, ">\n");  
       MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);  
       MPI_Wait(&req,&status);  
     }  
   }  
456    
457    // write actual data                  case Line3:
458    if(write_pointdata)                  case Tri6Face:
459    {                  case Rec8Face:
460      for (i_data =0 ;i_data<num_data;++i_data)                  case Line3_Contact:
461      {                  case Tri6Face_Contact:
462        if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])                  case Rec8Face_Contact:
463        {                      cellType = VTK_QUADRATIC_EDGE;
464          numPointsPerSample = elements->ReferenceElement->numQuadNodes;                      numVTKNodesPerElement = 3;
465          rank = getDataPointRank(data_pp[i_data]);                  break;
         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;  
         }  
466    
467          if( myRank == 0)                  case Tri6:
468          {          case Tri6Macro:    
469            char header[250];                  case Tet10Face:
470            sprintf(header, "<DataArray Name=\"%s\" type=\"Float32\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);                  case Tri6_Contact:
471            MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);                  case Tet10Face_Contact:
472            MPI_Wait(&req,&status);                      cellType = VTK_QUADRATIC_TRIANGLE;
473          }                      numVTKNodesPerElement = 6;
         // write out the data  
         // if the number of required components is more than the number  
         // of actual components, pad with zeros  
   
         char tmpbuf[15];  
         char* largebuf = MEMALLOC(nCompReqd*numLocalNodes*15 + numLocal*nCompReqd + nCompReqd + 15,char);  
         largebuf[0] = '\0';  
         bool_t do_write=TRUE;  
         size_t tsz = 0;  
   
         for(k=0;k < nodeCache.size;k++)  
         {  
           i = nodeCache.values[k];  
   
           if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)  
           {  
             if (mesh_p->Nodes->toReduced[i]>=0)  
             {  
               switch(getFunctionSpaceType(data_pp[i_data]))  
               {  
               case FINLEY_DEGREES_OF_FREEDOM:  
                 values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);  
474                  break;                  break;
475                case FINLEY_REDUCED_DEGREES_OF_FREEDOM:  
476                  values = getSampleData(data_pp[i_data],mesh_p->Nodes->reducedDegreeOfFreedom[i]);  
477                    case Rec8:
478                    case Hex20Face:
479                    case Rec8_Contact:
480                    case Hex20Face_Contact:
481                        cellType = VTK_QUADRATIC_QUAD;
482                        numVTKNodesPerElement = 8;
483                  break;                  break;
484                case FINLEY_NODES:  
485                  values = getSampleData(data_pp[i_data],i);          case Tet10Macro:
486                    case Tet10:
487                        cellType = VTK_QUADRATIC_TETRA;
488                        numVTKNodesPerElement = 10;
489                  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++)  
               {  
490    
491                  sprintf(tmpbuf," %e",values[m]);                  case Hex20:
492                  __STRCAT(largebuf,tmpbuf,tsz)                      cellType = VTK_QUADRATIC_HEXAHEDRON;
493                }                      numVTKNodesPerElement = 20;
494                for (m=0; m<nCompReqd-shape; m++)                  break;
495                {  
496                  __STRCAT(largebuf,zero,tsz)          case Hex27Macro:
497                }                  case Hex27:
498              }                      numCellFactor = 8;
499              else if (nCompReqd == 9)                      cellType = VTK_HEXAHEDRON;
500              {                      numVTKNodesPerElement = 8;
501                // tensor data, so have a 3x3 matrix to output as a row                  break;
502                //  of 9 data points  
503                count = 0;                  default:
504                for (m=0; m<shape; m++)                      sprintf(errorMsg, "saveVTK: Element type %s is not supported by VTK.", elements->referenceElementSet->referenceElement->Type->Name);
505                {                      Finley_setError(VALUE_ERROR, errorMsg);
506                  for (n=0; n<shape; n++)              }
507                  {          }
508                    sprintf(tmpbuf," %e", values[count]);      }
509                    __STRCAT(largebuf,tmpbuf,tsz)  
510                    count++;      /* allocate enough memory for text buffer */
511    
512        txtBufferSize = strlen(vtkHeader) + 3*LEN_INT_FORMAT + (30+3*maxNameLen)+strlen(metadata)+strlen(metadata_schema);
513        if (mpi_size > 1) {
514            txtBufferSize = MAX(txtBufferSize, myNumPoints * LEN_TMP_BUFFER);
515            txtBufferSize = MAX(txtBufferSize, numCellFactor * myNumCells *
516                    (LEN_INT_FORMAT * numVTKNodesPerElement + 1));
517            txtBufferSize = MAX(txtBufferSize,
518                    numCellFactor * myNumCells * LEN_TENSOR_FORMAT);
519            txtBufferSize = MAX(txtBufferSize, myNumPoints * LEN_TENSOR_FORMAT);
520        }
521        txtBuffer = TMPMEMALLOC(txtBufferSize+1, char);
522    
523        /* sets error if memory allocation failed */
524        Finley_checkPtr(txtBuffer);
525    
526        /************************************************************************/
527        /* write number of points and the mesh component */
528    
529        if (Finley_noError()) {
530            const index_t *nodeIndex;
531            if (FINLEY_REDUCED_NODES == nodeType) {
532                nodeIndex = elements->referenceElementSet->referenceElement->Type->linearNodes;
533            } else if (Line3Macro == typeId) {
534                 nodeIndex = VTK_LINE3_INDEX;
535            } else if ( (Rec9 == typeId) || (Rec9Macro == typeId) ) {
536                nodeIndex = VTK_REC9_INDEX;
537            } else if (Hex20 == typeId) {
538                nodeIndex = VTK_HEX20_INDEX;
539            } else if ( (Hex27 == typeId) || (Hex27Macro == typeId) ){
540                nodeIndex = VTK_HEX27_INDEX;
541            } else if (numVTKNodesPerElement  !=  elements->referenceElementSet->referenceElement->Type->numNodes) {
542                nodeIndex = elements->referenceElementSet->referenceElement->Type->relevantGeoNodes;
543            } else {
544                nodeIndex = NULL;
545            }
546    
547            if (strlen(metadata)>0) {
548               if (strlen(metadata_schema)>0) {
549                  sprintf(txtBuffer, vtkHeader," ",metadata_schema,metadata,"\n",globalNumPoints, numCellFactor*globalNumCells, 3);
550               } else {
551                  sprintf(txtBuffer, vtkHeader,"","",metadata,"\n",globalNumPoints, numCellFactor*globalNumCells, 3);
552               }
553            } else {
554               if (strlen(metadata_schema)>0) {
555                  sprintf(txtBuffer, vtkHeader," ",metadata_schema,"","",globalNumPoints, numCellFactor*globalNumCells, 3);
556               } else {
557                  sprintf(txtBuffer, vtkHeader,"","","","",globalNumPoints, numCellFactor*globalNumCells, 3);
558               }
559            }
560    
561            if (mpi_size > 1) {
562                /* write the nodes */
563                MPI_RANK0_WRITE_SHARED(txtBuffer);
564                txtBuffer[0] = '\0';
565                txtBufferInUse = 0;
566                if (nDim==2) {
567                    for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
568                        if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
569                            sprintf(tmpBuffer, VECTOR_FORMAT,
570                                mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
571                                mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
572                                0.);
573                            __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
574                        }
575                  }                  }
576                  for (n=0; n<3-shape; n++)              } else {
577                  {                  for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
578                    __STRCAT(largebuf,zero,tsz)                      if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
579                            sprintf(tmpBuffer, VECTOR_FORMAT,
580                                mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
581                                mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
582                                mesh_p->Nodes->Coordinates[INDEX2(2, i, nDim)]);
583                            __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
584                        }
585                  }                  }
586                }              } /* nDim */
587                for (m=0; m<3-shape; m++)              MPI_WRITE_ORDERED(txtBuffer);
588                {  
589                  for (n=0; n<3; n++)              /* write the cells */
590                  {              MPI_RANK0_WRITE_SHARED(tags_End_Points_and_Start_Conn);
591                    __STRCAT(largebuf,zero,tsz)              txtBuffer[0] = '\0';
592                txtBufferInUse = 0;
593                if (nodeIndex == NULL) {
594                    for (i = 0; i < numCells; i++) {
595                        if (elements->Owner[i] == my_mpi_rank) {
596                            for (j = 0; j < numVTKNodesPerElement; j++) {
597                                sprintf(tmpBuffer, INT_FORMAT, globalNodeIndex[elements->Nodes[INDEX2(j, i, NN)]]);
598                                __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
599                            }
600                            __STRCAT(txtBuffer, NEWLINE, txtBufferInUse);
601                        }
602                  }                  }
603                }              } else {
604              }                  for (i = 0; i < numCells; i++) {
605              __STRCAT(largebuf,newline,tsz)                      if (elements->Owner[i] == my_mpi_rank) {
606            }                          for (l = 0; l < numCellFactor; l++) {
607                                const int* idx=&nodeIndex[l*numVTKNodesPerElement];
608                                for (j = 0; j < numVTKNodesPerElement; j++) {
609                                    sprintf(tmpBuffer, INT_FORMAT, globalNodeIndex[elements->Nodes[INDEX2(idx[j], i, NN)]]);
610                                    __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
611                                }
612                                __STRCAT(txtBuffer, NEWLINE, txtBufferInUse);
613                            }
614                        }
615                    }
616                } /* nodeIndex */
617                MPI_WRITE_ORDERED(txtBuffer);
618    
619          }              /* write the offsets */
620          // Write out local data              MPI_RANK0_WRITE_SHARED(tags_End_Conn_and_Start_Offset);
621                txtBuffer[0] = '\0';
622                txtBufferInUse = 0;
623                for (i = numVTKNodesPerElement*(myFirstCell*numCellFactor+1);
624                        i <= (myFirstCell+myNumCells)*numVTKNodesPerElement*numCellFactor; i += numVTKNodesPerElement)
625                {
626                    sprintf(tmpBuffer, INT_NEWLINE_FORMAT, i);
627                    __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
628                }
629                MPI_WRITE_ORDERED(txtBuffer);
630    
631                /* write element type */
632                sprintf(tmpBuffer, INT_NEWLINE_FORMAT, cellType);
633                MPI_RANK0_WRITE_SHARED(tags_End_Offset_and_Start_Type);
634                txtBuffer[0] = '\0';
635                txtBufferInUse = 0;
636                for (i = numVTKNodesPerElement*(myFirstCell*numCellFactor+1);
637                        i <= (myFirstCell+myNumCells)*numVTKNodesPerElement*numCellFactor; i += numVTKNodesPerElement)
638                {
639                    __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
640                }
641                MPI_WRITE_ORDERED(txtBuffer);
642                /* finalize cell information */
643                strcpy(txtBuffer, "</DataArray>\n</Cells>\n");
644                MPI_RANK0_WRITE_SHARED(txtBuffer);
645    
646            } else { /***** mpi_size == 1 *****/
647    
648                /* write the nodes */
649                fputs(txtBuffer, fileHandle_p);
650                if (nDim==2) {
651                    for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
652                        if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
653                            fprintf(fileHandle_p, VECTOR_FORMAT,
654                                mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
655                                mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
656                                0.);
657                        }
658                    }
659                } else {
660                    for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
661                        if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
662                            fprintf(fileHandle_p, VECTOR_FORMAT,
663                                mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
664                                mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
665                                mesh_p->Nodes->Coordinates[INDEX2(2, i, nDim)]);
666                        }
667                    }
668                } /* nDim */
669    
670                /* write the cells */
671                fputs(tags_End_Points_and_Start_Conn, fileHandle_p);
672                if (nodeIndex == NULL) {
673                    for (i = 0; i < numCells; i++) {
674                        for (j = 0; j < numVTKNodesPerElement; j++) {
675                            fprintf(fileHandle_p, INT_FORMAT, globalNodeIndex[elements->Nodes[INDEX2(j, i, NN)]]);
676                        }
677                        fprintf(fileHandle_p, NEWLINE);
678                    }
679                } else {
680                    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                                fprintf(fileHandle_p, INT_FORMAT, globalNodeIndex[elements->Nodes[INDEX2(idx[j], i, NN)]]);
685                            }
686                            fprintf(fileHandle_p, NEWLINE);
687                        }
688                    }
689                } /* nodeIndex */
690    
691          largebuf[tsz] = '\0';              /* write the offsets */
692          MPI_File_write_ordered(fh,largebuf,tsz,MPI_CHAR,&status);              fputs(tags_End_Conn_and_Start_Offset, fileHandle_p);
693          MEMFREE(largebuf);              for (i = numVTKNodesPerElement; i <= numCells*numVTKNodesPerElement*numCellFactor; i += numVTKNodesPerElement) {
694          if( myRank == 0)                  fprintf(fileHandle_p, INT_NEWLINE_FORMAT, i);
695          {              }
696            char *tag = "</DataArray>\n";  
697            MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);              /* write element type */
698            MPI_Wait(&req,&status);              sprintf(tmpBuffer, INT_NEWLINE_FORMAT, cellType);
699                fputs(tags_End_Offset_and_Start_Type, fileHandle_p);
700                for (i = 0; i < numCells*numCellFactor; i++)
701                    fputs(tmpBuffer, fileHandle_p);
702                /* finalize cell information */
703                fputs("</DataArray>\n</Cells>\n", fileHandle_p);
704            } /* mpi_size */
705    
706        } /* Finley_noError */
707    
708        /************************************************************************/
709        /* write cell data */
710    
711        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                if (!Finley_noError())
751                    break;
752          }          }
       }  
753      }      }
754      // Finish off with closing tag      /* only continue if no error occurred */
755      if(myRank == 0)      if (writeCellData && Finley_noError()) {
756      {          strcat(txtBuffer, ">\n");
757        char* tag =  "</PointData>\n";          if ( mpi_size > 1) {
758        MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);              MPI_RANK0_WRITE_SHARED(txtBuffer);
759        MPI_Wait(&req,&status);          } else {
760      }              fputs(txtBuffer, fileHandle_p);
761      MPIO_DEBUG(" Done Writing Pointdata ")          }
762    }  
763    // end write_pointdata          /* write the arrays */
764            for (dataIdx = 0; dataIdx < num_data; dataIdx++) {
765    // tag and bag...                if (!isEmpty(data_pp[dataIdx]) && isCellCentered[dataIdx]) {
766    if (myRank == 0)                  dim_t numPointsPerSample=getNumDataPointsPerSample(data_pp[dataIdx]);
767    {                  dim_t rank = getDataPointRank(data_pp[dataIdx]);
768      char *footer = "</Piece>\n</UnstructuredGrid>\n</VTKFile>";                  dim_t nComp = getDataPointSize(data_pp[dataIdx]);
769      MPI_File_iwrite_shared(fh,footer,strlen(footer),MPI_CHAR,&req);                  dim_t nCompReqd = 1; /* number of components mpi_required */
770      MPI_Wait(&req,&status);                  if (rank == 0) {
771    }                      nCompReqd = 1;
772                        shape = 0;
773    MEMFREE(nodesGlobal);                  } else if (rank == 1) {
774    MEMFREE(nodeCache.values);                      nCompReqd = 3;
775    MEMFREE(elementCache.values);                      shape = getDataPointShape(data_pp[dataIdx], 0);
776  #ifdef MPIO_HINTS                      if (shape > 3) {
777    MPI_Info_free(&infoHints);                          Finley_setError(VALUE_ERROR, "saveVTK: rank 1 objects must have 3 components at most.");
778  #undef MPIO_HINTS                      }
779  #endif                  } else {
780    MPI_File_close(&fh);                      nCompReqd = 9;
781    MPIO_DEBUG(" ***** Exit saveVTK ***** ")                      shape = getDataPointShape(data_pp[dataIdx], 0);
782  #undef __STRCAT                      if (shape > 3 || shape != getDataPointShape(data_pp[dataIdx], 1)) {
783  }                          Finley_setError(VALUE_ERROR, "saveVTK: rank 2 objects of shape 2x2 or 3x3 supported only.");
784                        }
785  #undef MPIO_DEBUG                  }
786  #else                  /* bail out if an error occurred */
787                    if (!Finley_noError())
788                        break;
789    
790                    sprintf(txtBuffer, tag_Float_DataArray, names_p[dataIdx], nCompReqd);
791                    if ( mpi_size > 1) {
792                        MPI_RANK0_WRITE_SHARED(txtBuffer);
793                    } else {
794                        fputs(txtBuffer, fileHandle_p);
795                    }
796    
797                    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  void Finley_Mesh_saveVTK(const char * filename_p, Finley_Mesh *mesh_p, const dim_t num_data,char* *names_p, escriptDataC* *data_pp)          strcpy(txtBuffer, "</CellData>\n");
879  {          if ( mpi_size > 1) {
880    char error_msg[LenErrorMsg_MAX];              MPI_RANK0_WRITE_SHARED(txtBuffer);
881    /* if there is no mesh we just return */          } else {
882    if (mesh_p==NULL) return;              fputs(txtBuffer, fileHandle_p);
883            }
884    int i, j, k, numVTKNodesPerElement,i_data,m, count, n, rank,shape, numPoints, cellType, numCells,      } /* if noError && writeCellData */
885    nDim, numPointsPerSample, nComp, nCompReqd;  
886        /************************************************************************/
887    index_t j2;      /* write point data */
888    double* values, rtmp;  
889    char elemTypeStr[32];      if (writePointData && Finley_noError()) {
890            /* mark the active data arrays */
891    /* open the file and check handle */          bool_t set_scalar=FALSE, set_vector=FALSE, set_tensor=FALSE;
892            strcpy(txtBuffer, "<PointData");
893    FILE * fileHandle_p = fopen(filename_p, "w");          for (dataIdx = 0; dataIdx<num_data; dataIdx++) {
894    if (fileHandle_p==NULL)              if (!isEmpty(data_pp[dataIdx]) && !isCellCentered[dataIdx]) {
895    {                  switch (getDataPointRank(data_pp[dataIdx])) {
896      sprintf(error_msg, "saveVTK: File %s could not be opened for writing.", filename_p);                      case 0:
897      Finley_setError(IO_ERROR,error_msg);                          if (!set_scalar) {
898      return;                              strcat(txtBuffer, " Scalars=\"");
899    }                              strcat(txtBuffer, names_p[dataIdx]);
900    /* find the mesh type to be written */                              strcat(txtBuffer, "\"");
901    int nodetype=FINLEY_DEGREES_OF_FREEDOM;                              set_scalar = TRUE;
902    int elementtype=FINLEY_UNKNOWN;                          }
903    bool_t isCellCentered[num_data],write_celldata=FALSE,write_pointdata=FALSE;                      break;
904    for (i_data=0;i_data<num_data;++i_data)                      case 1:
905    {                          if (!set_vector) {
906      if (! isEmpty(data_pp[i_data]))                              strcat(txtBuffer, " Vectors=\"");
907      {                              strcat(txtBuffer, names_p[dataIdx]);
908        switch(getFunctionSpaceType(data_pp[i_data]))                              strcat(txtBuffer, "\"");
909        {                              set_vector = TRUE;
910        case FINLEY_DEGREES_OF_FREEDOM:                          }
911          nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;                      break;
912          if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)                      case 2:
913          {                          if (!set_tensor) {
914            elementtype=FINLEY_ELEMENTS;                              strcat(txtBuffer, " Tensors=\"");
915          }                              strcat(txtBuffer, names_p[dataIdx]);
916          else                              strcat(txtBuffer, "\"");
917          {                              set_tensor = TRUE;
918            Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");                          }
919            fclose(fileHandle_p);                      break;
920            return;                      default:
921          }                          sprintf(errorMsg, "saveVTK: data %s: VTK supports data with rank <= 2 only.", names_p[dataIdx]);
922          isCellCentered[i_data]=FALSE;                          Finley_setError(VALUE_ERROR, errorMsg);
923          break;                  }
924        case FINLEY_REDUCED_DEGREES_OF_FREEDOM:              }
925          nodetype = FINLEY_REDUCED_DEGREES_OF_FREEDOM;              if (!Finley_noError())
926          if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)                  break;
         {  
           elementtype=FINLEY_ELEMENTS;  
         }  
         else  
         {  
           Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");  
           fclose(fileHandle_p);  
           return;  
         }  
         isCellCentered[i_data]=FALSE;  
         break;  
       case FINLEY_NODES:  
         nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;  
         if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)  
         {  
           elementtype=FINLEY_ELEMENTS;  
         }  
         else  
         {  
           Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");  
           fclose(fileHandle_p);  
           return;  
         }  
         isCellCentered[i_data]=FALSE;  
         break;  
       case FINLEY_ELEMENTS:  
         nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;  
         if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)  
         {  
           elementtype=FINLEY_ELEMENTS;  
         }  
         else  
         {  
           Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");  
           fclose(fileHandle_p);  
           return;  
         }  
         isCellCentered[i_data]=TRUE;  
         break;  
       case FINLEY_FACE_ELEMENTS:  
         nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;  
         if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_FACE_ELEMENTS)  
         {  
           elementtype=FINLEY_FACE_ELEMENTS;  
         }  
         else  
         {  
           Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");  
           fclose(fileHandle_p);  
           return;  
         }  
         isCellCentered[i_data]=TRUE;  
         break;  
       case FINLEY_POINTS:  
         nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;  
         if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_POINTS)  
         {  
           elementtype=FINLEY_POINTS;  
         }  
         else  
         {  
           Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");  
           fclose(fileHandle_p);  
           return;  
         }  
         isCellCentered[i_data]=TRUE;  
         break;  
       case FINLEY_CONTACT_ELEMENTS_1:  
         nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;  
         if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1)  
         {  
           elementtype=FINLEY_CONTACT_ELEMENTS_1;  
         }  
         else  
         {  
           Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");  
           fclose(fileHandle_p);  
           return;  
         }  
         isCellCentered[i_data]=TRUE;  
         break;  
       case FINLEY_CONTACT_ELEMENTS_2:  
         nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;  
         if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1)  
         {  
           elementtype=FINLEY_CONTACT_ELEMENTS_1;  
         }  
         else  
         {  
           Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");  
           fclose(fileHandle_p);  
           return;  
927          }          }
         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=\"Float32\" 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");  
928      }      }
929        /* only continue if no error occurred */
930        if (writePointData && Finley_noError()) {
931            strcat(txtBuffer, ">\n");
932            if ( mpi_size > 1) {
933                MPI_RANK0_WRITE_SHARED(txtBuffer);
934            } else {
935                fputs(txtBuffer, fileHandle_p);
936            }
937    
938            /* write the arrays */
939            for (dataIdx = 0; dataIdx < num_data; dataIdx++) {
940                if (!isEmpty(data_pp[dataIdx]) && !isCellCentered[dataIdx]) {
941                    Finley_NodeMapping* nodeMapping;
942                    dim_t rank = getDataPointRank(data_pp[dataIdx]);
943                    dim_t nCompReqd = 1; /* number of components mpi_required */
944                    if (getFunctionSpaceType(data_pp[dataIdx]) == FINLEY_REDUCED_NODES) {
945                        nodeMapping = mesh_p->Nodes->reducedNodesMapping;
946                    } else {
947                        nodeMapping = mesh_p->Nodes->nodesMapping;
948                    }
949                    if (rank == 0) {
950                        nCompReqd = 1;
951                        shape = 0;
952                    } else if (rank == 1) {
953                        nCompReqd = 3;
954                        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    fprintf(fileHandle_p, "</DataArray>\n");                  txtBufferInUse = 0;
978    fprintf(fileHandle_p, "</Points>\n");                  for (i=0; i<mesh_p->Nodes->numNodes; i++) {
979                        k = globalNodeIndex[i];
980    /* write out the DataArray element for the connectivity */                      if ( (myFirstNode <= k) && (k < myLastNode) ) {
981                            __const double *values = getSampleDataRO(data_pp[dataIdx], nodeMapping->target[i]);
982    int NN = elements->ReferenceElement->Type->numNodes;                          /* if the number of mpi_required components is more than
983    fprintf(fileHandle_p, "<Cells>\n");                           * the number of actual components, pad with zeros.
984    fprintf(fileHandle_p, "<DataArray Name=\"connectivity\" type=\"Int32\" format=\"ascii\">\n");                           * Probably only need to get shape of first element */
985                            if (nCompReqd == 1) {
986    if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)                              sprintf(tmpBuffer, SCALAR_FORMAT, values[0]);
987    {                          } else if (nCompReqd == 3) {
988      for (i = 0; i < numCells; i++)                              if (shape==1) {
989      {                                  sprintf(tmpBuffer, VECTOR_FORMAT,
990        for (j = 0; j < numVTKNodesPerElement; j++)                                          values[0], 0.f, 0.f);
991          fprintf(fileHandle_p,"%d ",mesh_p->Nodes->toReduced[elements->Nodes[INDEX2(elements->ReferenceElement->Type->linearNodes[j], i, NN)]]);                              } else if (shape==2) {
992        fprintf(fileHandle_p, "\n");                                  sprintf(tmpBuffer, VECTOR_FORMAT,
993      }                                          values[0], values[1], 0.f);
994    }                              } else if (shape==3) {
995    else if (VTK_QUADRATIC_HEXAHEDRON==cellType)                                  sprintf(tmpBuffer, VECTOR_FORMAT,
996    {                                          values[0], values[1], values[2]);
997      for (i = 0; i < numCells; i++)                              }
998      {                          } else if (nCompReqd == 9) {
999        fprintf(fileHandle_p,"%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d\n",                              if (shape==1) {
1000                elements->Nodes[INDEX2(0, i, NN)],                                  sprintf(tmpBuffer, TENSOR_FORMAT,
1001                elements->Nodes[INDEX2(1, i, NN)],                                      values[0], 0.f, 0.f,
1002                elements->Nodes[INDEX2(2, i, NN)],                                            0.f, 0.f, 0.f,
1003                elements->Nodes[INDEX2(3, i, NN)],                                            0.f, 0.f, 0.f);
1004                elements->Nodes[INDEX2(4, i, NN)],                              } else if (shape==2) {
1005                elements->Nodes[INDEX2(5, i, NN)],                                  sprintf(tmpBuffer, TENSOR_FORMAT,
1006                elements->Nodes[INDEX2(6, i, NN)],                                      values[0], values[1], 0.f,
1007                elements->Nodes[INDEX2(7, i, NN)],                                      values[2], values[3], 0.f,
1008                elements->Nodes[INDEX2(8, i, NN)],                                            0.f,       0.f, 0.f);
1009                elements->Nodes[INDEX2(9, i, NN)],                              } else if (shape==3) {
1010                elements->Nodes[INDEX2(10, i, NN)],                                  sprintf(tmpBuffer, TENSOR_FORMAT,
1011                elements->Nodes[INDEX2(11, i, NN)],                                      values[0], values[1], values[2],
1012                elements->Nodes[INDEX2(16, i, NN)],                                      values[3], values[4], values[5],
1013                elements->Nodes[INDEX2(17, i, NN)],                                      values[6], values[7], values[8]);
1014                elements->Nodes[INDEX2(18, i, NN)],                              }
1015                elements->Nodes[INDEX2(19, i, NN)],                          }
1016                elements->Nodes[INDEX2(12, i, NN)],                          if ( mpi_size > 1) {
1017                elements->Nodes[INDEX2(13, i, NN)],                              __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
1018                elements->Nodes[INDEX2(14, i, NN)],                          } else {
1019                elements->Nodes[INDEX2(15, i, NN)]);                              fputs(tmpBuffer, fileHandle_p);
1020      }                          }
1021    }                      } /* if this is my node */
1022    else if (numVTKNodesPerElement!=NN)                  } /* for i (numNodes) */
1023    {  
1024      for (i = 0; i < numCells; i++)                  if ( mpi_size > 1) {
1025      {                      MPI_WRITE_ORDERED(txtBuffer);
1026        for (j = 0; j < numVTKNodesPerElement; j++) fprintf(fileHandle_p,"%d ", elements->Nodes[INDEX2(elements->ReferenceElement->Type->geoNodes[j], i, NN)]);                      MPI_RANK0_WRITE_SHARED(tag_End_DataArray);
1027        fprintf(fileHandle_p, "\n");                  } else {
1028      }                      fputs(tag_End_DataArray, fileHandle_p);
1029    }                  }
1030    else              } /* !isEmpty && !isCellCentered */
1031    {          } /* for dataIdx */
     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");  
1032    
1033    /* write out the DataArray element for the offsets */          strcpy(txtBuffer, "</PointData>\n");
1034    fprintf(fileHandle_p, "<DataArray Name=\"offsets\" type=\"Int32\" format=\"ascii\">\n");          if ( mpi_size > 1) {
1035    for (i=numVTKNodesPerElement; i<=numCells*numVTKNodesPerElement; i+=numVTKNodesPerElement) fprintf(fileHandle_p, "%d\n", i);              MPI_RANK0_WRITE_SHARED(txtBuffer);
1036    fprintf(fileHandle_p, "</DataArray>\n");          } else {
1037                fputs(txtBuffer, fileHandle_p);
   /* write out the DataArray element for the types */  
   fprintf(fileHandle_p, "<DataArray Name=\"types\" type=\"UInt8\" format=\"ascii\">\n");  
   for (i=0; i<numCells; i++) fprintf(fileHandle_p, "%d\n", cellType);  
   fprintf(fileHandle_p, "</DataArray>\n");  
   
   /* finish off the <Cells> element */  
   fprintf(fileHandle_p, "</Cells>\n");  
   
   /* cell data */  
   if (write_celldata)  
   {  
     /* mark the active data arrays */  
     bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;  
     fprintf(fileHandle_p, "<CellData");  
     for (i_data =0 ;i_data<num_data;++i_data)  
     {  
       if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data])  
       {  
         /* if the rank == 0:   --> scalar data  
         * if the rank == 1:   --> vector data  
         * if the rank == 2:   --> tensor data  
         */  
         switch(getDataPointRank(data_pp[i_data]))  
         {  
         case 0:  
           if (! set_scalar)  
           {  
             fprintf(fileHandle_p," Scalars=\"%s\"",names_p[i_data]);  
             set_scalar=TRUE;  
           }  
           break;  
         case 1:  
           if (! set_vector)  
           {  
             fprintf(fileHandle_p," Vectors=\"%s\"",names_p[i_data]);  
             set_vector=TRUE;  
           }  
           break;  
         case 2:  
           if (! set_tensor)  
           {  
             fprintf(fileHandle_p," Tensors=\"%s\"",names_p[i_data]);  
             set_tensor=TRUE;  
           }  
           break;  
         default:  
           sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);  
           Finley_setError(VALUE_ERROR,error_msg);  
           fclose(fileHandle_p);  
           return;  
1038          }          }
1039        }      } /* if noError && writePointData */
     }  
     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=\"Float32\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);  
1040    
1041          double sampleAvg[nComp];      /* Final write to VTK file */
1042          for (i=0; i<numCells; i++)      if (Finley_noError()) {
1043          {          if (mpi_size > 1) {
1044            values = getSampleData(data_pp[i_data], i);              MPI_RANK0_WRITE_SHARED(vtkFooter);
1045            /* averaging over the number of points in the sample */          } else {
1046            for (k=0; k<nComp; k++)              fputs(vtkFooter, fileHandle_p);
           {  
             rtmp = 0.;  
             for (j=0; j<numPointsPerSample; j++) rtmp += values[INDEX2(k,j,nComp)];  
             sampleAvg[k] = rtmp/numPointsPerSample;  
           }  
           /* if the number of required components is more than the number  
           * of actual components, pad with zeros  
           */  
           /* probably only need to get shape of first element */  
           /* write the data different ways for scalar, vector and tensor */  
           if (nCompReqd == 1)  
           {  
             fprintf(fileHandle_p, " %e", sampleAvg[0]);  
           }  
           else if (nCompReqd == 3)  
           {  
             /* write out the data */  
             for (m=0; m<shape; m++) fprintf(fileHandle_p, " %e", sampleAvg[m]);  
             for (m=0; m<nCompReqd-shape; m++) fprintf(fileHandle_p, " %e", 0.);  
           }  
           else if (nCompReqd == 9)  
           {  
             /* tensor data, so have a 3x3 matrix to output as a row  
             * of 9 data points */  
             count = 0;  
             for (m=0; m<shape; m++)  
             {  
               for (n=0; n<shape; n++)  
               {  
                 fprintf(fileHandle_p, " %e", sampleAvg[count]);  
                 count++;  
               }  
               for (n=0; n<3-shape; n++) fprintf(fileHandle_p, " %e", 0.);  
             }  
             for (m=0; m<3-shape; m++)  
               for (n=0; n<3; n++) fprintf(fileHandle_p, " %e", 0.);  
           }  
           fprintf(fileHandle_p, "\n");  
1047          }          }
         fprintf(fileHandle_p, "</DataArray>\n");  
       }  
1048      }      }
1049      fprintf(fileHandle_p, "</CellData>\n");  
1050    }      if ( mpi_size > 1) {
1051    /* point data */  #ifdef PASO_MPI
1052    if (write_pointdata)          MPI_File_close(&mpi_fileHandle_p);
1053    {          MPI_Barrier(mesh_p->Nodes->MPIInfo->comm);
1054      /* mark the active data arrays */  #endif
1055      bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;      } else {
1056      fprintf(fileHandle_p, "<PointData");          fclose(fileHandle_p);
     for (i_data =0 ;i_data<num_data;++i_data)  
     {  
       if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])  
       {  
         /* if the rank == 0:   --> scalar data  
         * if the rank == 1:   --> vector data  
         * if the rank == 2:   --> tensor data  
         */  
         switch(getDataPointRank(data_pp[i_data]))  
         {  
         case 0:  
           if (! set_scalar)  
           {  
             fprintf(fileHandle_p," Scalars=\"%s\"",names_p[i_data]);  
             set_scalar=TRUE;  
           }  
           break;  
         case 1:  
           if (! set_vector)  
           {  
             fprintf(fileHandle_p," Vectors=\"%s\"",names_p[i_data]);  
             set_vector=TRUE;  
           }  
           break;  
         case 2:  
           if (! set_tensor)  
           {  
             fprintf(fileHandle_p," Tensors=\"%s\"",names_p[i_data]);  
             set_tensor=TRUE;  
           }  
           break;  
         default:  
           sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);  
           Finley_setError(VALUE_ERROR,error_msg);  
           fclose(fileHandle_p);  
           return;  
         }  
       }  
     }  
     fprintf(fileHandle_p, ">\n");  
     /* write the arrays */  
     for (i_data =0 ;i_data<num_data;++i_data)  
     {  
       if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])  
       {  
         numPointsPerSample = elements->ReferenceElement->numQuadNodes;  
         rank = getDataPointRank(data_pp[i_data]);  
         nComp = getDataPointSize(data_pp[i_data]);  
         nCompReqd=1;   /* the number of components required by vtk */  
         shape=0;  
         if (rank == 0)  
         {  
           nCompReqd = 1;  
         }  
         else if (rank == 1)  
         {  
           shape=getDataPointShape(data_pp[i_data], 0);  
           if  (shape>3)  
           {  
             Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");  
             fclose(fileHandle_p);  
             return;  
           }  
           nCompReqd = 3;  
         }  
         else  
         {  
           shape=getDataPointShape(data_pp[i_data], 0);  
           if  (shape>3 || shape != getDataPointShape(data_pp[i_data], 1))  
           {  
             Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");  
             fclose(fileHandle_p);  
             return;  
           }  
           nCompReqd = 9;  
         }  
         fprintf(fileHandle_p, "<DataArray Name=\"%s\" type=\"Float32\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);  
         /* write out the data */  
         /* if the number of required components is more than the number  
         * of actual components, pad with zeros  
         */  
         bool_t do_write=TRUE;  
         for (i=0; i<mesh_p->Nodes->numNodes; i++)  
         {  
           if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)  
           {  
             if (mesh_p->Nodes->toReduced[i]>=0)  
             {  
               switch(getFunctionSpaceType(data_pp[i_data]))  
               {  
               case FINLEY_DEGREES_OF_FREEDOM:  
                 values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);  
                 break;  
               case FINLEY_REDUCED_DEGREES_OF_FREEDOM:  
                 values = getSampleData(data_pp[i_data],mesh_p->Nodes->reducedDegreeOfFreedom[i]);  
                 break;  
               case FINLEY_NODES:  
                 values = getSampleData(data_pp[i_data],i);  
                 break;  
               }  
               do_write=TRUE;  
             }  
             else  
             {  
               do_write=FALSE;  
             }  
           }  
           else  
           {  
             do_write=TRUE;  
             switch(getFunctionSpaceType(data_pp[i_data]))  
             {  
             case FINLEY_DEGREES_OF_FREEDOM:  
               values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);  
               break;  
             case FINLEY_NODES:  
               values = getSampleData(data_pp[i_data],i);  
               break;  
             }  
           }  
           /* write the data different ways for scalar, vector and tensor */  
           if (do_write)  
           {  
             if (nCompReqd == 1)  
             {  
               fprintf(fileHandle_p, " %e", values[0]);  
             }  
             else if (nCompReqd == 3)  
             {  
               for (m=0; m<shape; m++) fprintf(fileHandle_p, " %e", values[m]);  
               for (m=0; m<nCompReqd-shape; m++) fprintf(fileHandle_p, " %e", 0.);  
             }  
             else if (nCompReqd == 9)  
             {  
               /* tensor data, so have a 3x3 matrix to output as a row  
               * of 9 data points */  
               count = 0;  
               for (m=0; m<shape; m++)  
               {  
                 for (n=0; n<shape; n++)  
                 {  
                   fprintf(fileHandle_p, " %e", values[count]);  
                   count++;  
                 }  
                 for (n=0; n<3-shape; n++) fprintf(fileHandle_p, " %e", 0.);  
               }  
               for (m=0; m<3-shape; m++)  
                 for (n=0; n<3; n++) fprintf(fileHandle_p, " %e", 0.);  
             }  
             fprintf(fileHandle_p, "\n");  
           }  
         }  
         fprintf(fileHandle_p, "</DataArray>\n");  
       }  
1057      }      }
1058      fprintf(fileHandle_p, "</PointData>\n");      TMPMEMFREE(isCellCentered);
1059    }      TMPMEMFREE(txtBuffer);
   /* finish off the piece */  
   fprintf(fileHandle_p, "</Piece>\n");  
   
   fprintf(fileHandle_p, "</UnstructuredGrid>\n");  
   /* write the xml footer */  
   fprintf(fileHandle_p, "</VTKFile>\n");  
   /* close the file */  
   fclose(fileHandle_p);  
   return;  
1060  }  }
 #endif  
1061    

Legend:
Removed from v.818  
changed lines
  Added in v.2814

  ViewVC Help
Powered by ViewVC 1.1.26