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

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

  ViewVC Help
Powered by ViewVC 1.1.26