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

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

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.26