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

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

  ViewVC Help
Powered by ViewVC 1.1.26