/[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

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

Legend:
Removed from v.794  
changed lines
  Added in v.1387

  ViewVC Help
Powered by ViewVC 1.1.26