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

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

  ViewVC Help
Powered by ViewVC 1.1.26