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

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

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.26