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

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

  ViewVC Help
Powered by ViewVC 1.1.26