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

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

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.26