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

Legend:
Removed from v.793  
changed lines
  Added in v.2141

  ViewVC Help
Powered by ViewVC 1.1.26