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

Legend:
Removed from v.794  
changed lines
  Added in v.2548

  ViewVC Help
Powered by ViewVC 1.1.26