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

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

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.26