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

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

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

trunk/finley/src/finley/Mesh_saveVTK.c revision 201 by jgs, Wed Nov 23 04:10:21 2005 UTC trunk/finley/src/Mesh_saveVTK.c revision 2141 by caltinay, Tue Dec 9 04:32:32 2008 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-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  /*   Author: Paul Cochrane, cochrane@esscc.uq.edu.au */  #else
 /*   Version: $Id$ */  
77    
78  /**************************************************************/  #define MPI_WRITE_ORDERED(A, B)
79    #define MPI_RANK0_WRITE_SHARED(A)
80    
81  #include "Mesh.h"  #endif /* PASO_MPI */
82  #include "vtkCellType.h"  /* copied from vtk source directory !!! */  
83    
84    /* Returns one if the node given by coords and idx is within the quadrant
85     * 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;
360                   elementtype=FINLEY_ELEMENTS;          }
361               } else {          if (elements==NULL) {
362                   Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");              Finley_setError(SYSTEM_ERROR, "saveVTK: undefined element file");
363                   fclose(fileHandle_p);          } else {
364                   return;              /* map finley element type to VTK element type */
365               }              numCells = elements->numElements;
366               isCellCentered[i_data]=FALSE;              globalNumCells = Finley_ElementFile_getGlobalNumElements(elements);
367               break;              myNumCells = Finley_ElementFile_getMyNumElements(elements);
368             case FINLEY_REDUCED_DEGREES_OF_FREEDOM:              myFirstCell = Finley_ElementFile_getFirstElement(elements);
369               nodetype = FINLEY_REDUCED_DEGREES_OF_FREEDOM;              NN = elements->numNodes;
370               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {              if (nodeType==FINLEY_REDUCED_NODES) {
371                   elementtype=FINLEY_ELEMENTS;                  typeId = elements->LinearReferenceElement->Type->TypeId;
372               } else {                  if (hasReducedElements) {
373                   Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");                      quadNodes_p=elements->LinearReferenceElementReducedOrder->QuadNodes;
374                   fclose(fileHandle_p);                  } else {
375                   return;                      quadNodes_p=elements->LinearReferenceElement->QuadNodes;
376               }                  }
377               isCellCentered[i_data]=FALSE;              } else {
378               break;                  typeId = elements->ReferenceElement->Type->TypeId;
379             case FINLEY_NODES:                  if (hasReducedElements)
380               nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;                      quadNodes_p=elements->ReferenceElementReducedOrder->QuadNodes;
381               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {                  else
382                   elementtype=FINLEY_ELEMENTS;                      quadNodes_p=elements->ReferenceElement->QuadNodes;
383               } else {              }
384                   Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");              switch (typeId) {
385                   fclose(fileHandle_p);                  case Point1:
386                   return;                  case Line2Face:
387               }                  case Line3Face:
388               isCellCentered[i_data]=FALSE;                  case Point1_Contact:
389               break;                  case Line2Face_Contact:
390             case FINLEY_ELEMENTS:                  case Line3Face_Contact:
391               nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;                      cellType = VTK_VERTEX;
392               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {                      numVTKNodesPerElement = 1;
393                   elementtype=FINLEY_ELEMENTS;                  break;
394               } else {  
395                   Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");                  case Line2:
396                   fclose(fileHandle_p);                  case Tri3Face:
397                   return;                  case Rec4Face:
398               }                  case Line2_Contact:
399               isCellCentered[i_data]=TRUE;                  case Tri3_Contact:
400               break;                  case Tri3Face_Contact:
401             case FINLEY_FACE_ELEMENTS:                  case Rec4Face_Contact:
402               nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;                      cellType = VTK_LINE;
403               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_FACE_ELEMENTS) {                      numVTKNodesPerElement = 2;
404                   elementtype=FINLEY_FACE_ELEMENTS;                  break;
405               } else {  
406                   Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");                  case Tri3:
407                   fclose(fileHandle_p);                  case Tet4Face:
408                   return;                  case Tet4Face_Contact:
409               }                      cellType = VTK_TRIANGLE;
410               isCellCentered[i_data]=TRUE;                      numVTKNodesPerElement = 3;
411               break;                  break;
412             case FINLEY_POINTS:  
413               nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;                  case Rec4:
414               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_POINTS) {                  case Hex8Face:
415                   elementtype=FINLEY_POINTS;                  case Rec4_Contact:
416               } else {                  case Hex8Face_Contact:
417                   Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");                      cellType = VTK_QUAD;
418                   fclose(fileHandle_p);                      numVTKNodesPerElement = 4;
419                   return;                  break;
420               }  
421               isCellCentered[i_data]=TRUE;                  case Rec9:
422               break;                      numCellFactor = 4;
423             case FINLEY_CONTACT_ELEMENTS_1:                      cellType = VTK_QUAD;
424               nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;                      numVTKNodesPerElement = 4;
425               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {                  break;
426                   elementtype=FINLEY_CONTACT_ELEMENTS_1;  
427               } else {                  case Tet4:
428                   Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");                      cellType = VTK_TETRA;
429                   fclose(fileHandle_p);                      numVTKNodesPerElement = 4;
430                   return;                  break;
431               }  
432               isCellCentered[i_data]=TRUE;                  case Hex8:
433               break;                      cellType = VTK_HEXAHEDRON;
434             case FINLEY_CONTACT_ELEMENTS_2:                      numVTKNodesPerElement = 8;
435               nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;                  break;
436               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {  
437                   elementtype=FINLEY_CONTACT_ELEMENTS_1;                  case Line3:
438               } else {                  case Tri6Face:
439                   Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");                  case Rec8Face:
440                   fclose(fileHandle_p);                  case Line3_Contact:
441                   return;                  case Tri6Face_Contact:
442               }                  case Rec8Face_Contact:
443               isCellCentered[i_data]=TRUE;                      cellType = VTK_QUADRATIC_EDGE;
444               break;                      numVTKNodesPerElement = 3;
445             default:                  break;
446               sprintf(error_msg,"saveVTK: Finley does not know anything about function space type %d",getFunctionSpaceType(data_pp[i_data]));  
447               Finley_setError(TYPE_ERROR,error_msg);                  case Tri6:
448               fclose(fileHandle_p);                  case Tet10Face:
449               return;                  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          if (isCellCentered[i_data]) {      }
485             write_celldata=TRUE;  
486        /* allocate enough memory for text buffer */
487    
488        txtBufferSize = strlen(vtkHeader) + 3*LEN_INT_FORMAT + (30+3*maxNameLen);
489    
490        if (mpi_size > 1) {
491            txtBufferSize = MAX(txtBufferSize, myNumPoints * LEN_TMP_BUFFER);
492            txtBufferSize = MAX(txtBufferSize, numCellFactor * myNumCells *
493                    (LEN_INT_FORMAT * numVTKNodesPerElement + 1));
494            txtBufferSize = MAX(txtBufferSize,
495                    numCellFactor * myNumCells * LEN_TENSOR_FORMAT);
496            txtBufferSize = MAX(txtBufferSize, myNumPoints * LEN_TENSOR_FORMAT);
497        }
498        txtBuffer = TMPMEMALLOC(txtBufferSize+1, char);
499    
500        /* sets error if memory allocation failed */
501        Finley_checkPtr(txtBuffer);
502    
503        /************************************************************************/
504        /* write number of points and the mesh component */
505    
506        if (Finley_noError()) {
507            const index_t *nodeIndex;
508            if (FINLEY_REDUCED_NODES == nodeType) {
509                nodeIndex = elements->ReferenceElement->Type->linearNodes;
510            } else if (Rec9 == typeId) {
511                nodeIndex = VTK_REC9_INDEX;
512            } else if (Hex20 == typeId) {
513                nodeIndex = VTK_HEX20_INDEX;
514            } else if (Hex27 == typeId) {
515                nodeIndex = VTK_HEX27_INDEX;
516            } else if (numVTKNodesPerElement != NN) {
517                nodeIndex = elements->ReferenceElement->Type->geoNodes;
518          } else {          } else {
519             write_pointdata=TRUE;              nodeIndex = NULL;
520          }          }
521       }  
522    }          sprintf(txtBuffer, vtkHeader, globalNumPoints,
523    /* select nomber of points and the mesh component */                  numCellFactor*globalNumCells, 3);
524    int numPoints = mesh_p->Nodes->numNodes;  
525    if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {          if (mpi_size > 1) {
526        numPoints = mesh_p->Nodes->reducedNumNodes;              /* write the nodes */
527    } else {              MPI_RANK0_WRITE_SHARED(txtBuffer);
528        numPoints = mesh_p->Nodes->numNodes;              txtBuffer[0] = '\0';
529    }              txtBufferInUse = 0;
530    if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;              if (nDim==2) {
531    Finley_ElementFile* elements=NULL;                  for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
532    switch(elementtype) {                      if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
533      case FINLEY_ELEMENTS:                          sprintf(tmpBuffer, VECTOR_FORMAT,
534        elements=mesh_p->Elements;                              mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
535        break;                              mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
536      case FINLEY_FACE_ELEMENTS:                              0.);
537        elements=mesh_p->FaceElements;                          __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
538        break;                      }
539      case FINLEY_POINTS:                  }
540        elements=mesh_p->Points;              } else {
541        break;                  for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
542      case FINLEY_CONTACT_ELEMENTS_1:                      if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
543        elements=mesh_p->ContactElements;                          sprintf(tmpBuffer, VECTOR_FORMAT,
544        break;                              mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
545    }                              mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
546    if (elements==NULL) {                              mesh_p->Nodes->Coordinates[INDEX2(2, i, nDim)]);
547       Finley_setError(SYSTEM_ERROR,"saveVTK: undefined element file");                          __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
548       fclose(fileHandle_p);                      }
549       return;                  }
550    }              } /* nDim */
551    /* map finley element type to VTK element type */              MPI_WRITE_ORDERED(txtBuffer, txtBufferInUse);
552    int numCells = elements->numElements;    
553    int cellType;              /* write the cells */
554    ElementTypeId TypeId;              MPI_RANK0_WRITE_SHARED(tags_End_Points_and_Start_Conn);
555    char elemTypeStr[32];              txtBuffer[0] = '\0';
556    if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {              txtBufferInUse = 0;
557       TypeId = elements->LinearReferenceElement->Type->TypeId;              if (nodeIndex == NULL) {
558    } else {                  for (i = 0; i < numCells; i++) {
559       TypeId = elements->ReferenceElement->Type->TypeId;                      if (elements->Owner[i] == my_mpi_rank) {
560    }                          for (j = 0; j < numVTKNodesPerElement; j++) {
561                                sprintf(tmpBuffer, INT_FORMAT, globalNodeIndex[elements->Nodes[INDEX2(j, i, NN)]]);
562    switch(TypeId) {                              __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
563      case Point1:                          }
564      case Line2Face:                          __STRCAT(txtBuffer, NEWLINE, txtBufferInUse);
565      case Line3Face:                      }
566      case Point1_Contact:                  }
567      case Line2Face_Contact:              } else {
568      case Line3Face_Contact:                  for (i = 0; i < numCells; i++) {
569        cellType = VTK_VERTEX;                      if (elements->Owner[i] == my_mpi_rank) {
570        numVTKNodesPerElement = 1;                          for (l = 0; l < numCellFactor; l++) {
571        strcpy(elemTypeStr, "VTK_VERTEX");                              const int* idx=&nodeIndex[l*numVTKNodesPerElement];
572        break;                              for (j = 0; j < numVTKNodesPerElement; j++) {
573                                    sprintf(tmpBuffer, INT_FORMAT, globalNodeIndex[elements->Nodes[INDEX2(idx[j], i, NN)]]);
574      case Line2:                                  __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
575      case Tri3Face:                              }
576      case Rec4Face:                              __STRCAT(txtBuffer, NEWLINE, txtBufferInUse);
577      case Line2_Contact:                          }
578      case Tri3_Contact:                      }
579      case Tri3Face_Contact:                  }
580      case Rec4Face_Contact:              } /* nodeIndex */
581        cellType = VTK_LINE;              MPI_WRITE_ORDERED(txtBuffer, txtBufferInUse);
582        numVTKNodesPerElement = 2;  
583        strcpy(elemTypeStr, "VTK_LINE");              /* write the offsets */
584        break;              MPI_RANK0_WRITE_SHARED(tags_End_Conn_and_Start_Offset);
585                txtBuffer[0] = '\0';
586      case Tri3:              txtBufferInUse = 0;
587      case Tet4Face:              for (i = numVTKNodesPerElement*(myFirstCell*numCellFactor+1);
588      case Tet4Face_Contact:                      i <= (myFirstCell+myNumCells)*numVTKNodesPerElement*numCellFactor; i += numVTKNodesPerElement)
589        cellType = VTK_TRIANGLE;              {
590        numVTKNodesPerElement = 3;                  sprintf(tmpBuffer, INT_NEWLINE_FORMAT, i);
591        strcpy(elemTypeStr, "VTK_TRIANGLE");                  __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
592        break;              }
593                MPI_WRITE_ORDERED(txtBuffer, txtBufferInUse);
594      case Rec4:  
595      case Hex8Face:              /* write element type */
596      case Rec4_Contact:              sprintf(tmpBuffer, INT_NEWLINE_FORMAT, cellType);
597      case Hex8Face_Contact:              MPI_RANK0_WRITE_SHARED(tags_End_Offset_and_Start_Type);
598        cellType = VTK_QUAD;              txtBuffer[0] = '\0';
599        numVTKNodesPerElement = 4;              txtBufferInUse = 0;
600        strcpy(elemTypeStr, "VTK_QUAD");              for (i = numVTKNodesPerElement*(myFirstCell*numCellFactor+1);
601        break;                      i <= (myFirstCell+myNumCells)*numVTKNodesPerElement*numCellFactor; i += numVTKNodesPerElement)
602                {
603      case Tet4:                  __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
604        cellType = VTK_TETRA;              }
605        numVTKNodesPerElement = 4;              MPI_WRITE_ORDERED(txtBuffer, txtBufferInUse);
606        strcpy(elemTypeStr, "VTK_TETRA");              /* finalize cell information */
607        break;              strcpy(txtBuffer, "</DataArray>\n</Cells>\n");
608                MPI_RANK0_WRITE_SHARED(txtBuffer);
609      case Hex8:  
610        cellType = VTK_HEXAHEDRON;          } else { /***** mpi_size == 1 *****/
611        numVTKNodesPerElement = 8;  
612        strcpy(elemTypeStr, "VTK_HEXAHEDRON");              /* write the nodes */
613        break;              fputs(txtBuffer, fileHandle_p);
614                if (nDim==2) {
615      case Line3:                  for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
616      case Tri6Face:                      if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
617      case Rec8Face:                          fprintf(fileHandle_p, VECTOR_FORMAT,
618      case Line3_Contact:                              mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
619      case Tri6Face_Contact:                              mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
620      case Rec8Face_Contact:                              0.);
621        cellType = VTK_QUADRATIC_EDGE;                      }
622        numVTKNodesPerElement = 3;                  }
623        strcpy(elemTypeStr, "VTK_QUADRATIC_EDGE");              } else {
624        break;                  for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
625                        if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
626      case Tri6:                          fprintf(fileHandle_p, VECTOR_FORMAT,
627      case Tet10Face:                              mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
628      case Tri6_Contact:                              mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
629      case Tet10Face_Contact:                              mesh_p->Nodes->Coordinates[INDEX2(2, i, nDim)]);
630        cellType = VTK_QUADRATIC_TRIANGLE;                      }
631        numVTKNodesPerElement = 6;                  }
632        strcpy(elemTypeStr, "VTK_QUADRATIC_TRIANGLE");              } /* nDim */
633        break;  
634                /* write the cells */
635      case Rec8:              fputs(tags_End_Points_and_Start_Conn, fileHandle_p);
636      case Hex20Face:              if (nodeIndex == NULL) {
637      case Rec8_Contact:                  for (i = 0; i < numCells; i++) {
638      case Hex20Face_Contact:                      for (j = 0; j < numVTKNodesPerElement; j++) {
639        cellType = VTK_QUADRATIC_QUAD;                          fprintf(fileHandle_p, INT_FORMAT, globalNodeIndex[elements->Nodes[INDEX2(j, i, NN)]]);
640        numVTKNodesPerElement = 8;                      }
641        strcpy(elemTypeStr, "VTK_QUADRATIC_QUAD");                      fprintf(fileHandle_p, NEWLINE);
642        break;                  }
643                } else {
644      case Tet10:                  for (i = 0; i < numCells; i++) {
645        cellType = VTK_QUADRATIC_TETRA;                      for (l = 0; l < numCellFactor; l++) {
646        numVTKNodesPerElement = 10;                          const int* idx=&nodeIndex[l*numVTKNodesPerElement];
647        strcpy(elemTypeStr, "VTK_QUADRATIC_TETRA");                          for (j = 0; j < numVTKNodesPerElement; j++) {
648        break;                              fprintf(fileHandle_p, INT_FORMAT, globalNodeIndex[elements->Nodes[INDEX2(idx[j], i, NN)]]);
649                            }
650      case Hex20:                          fprintf(fileHandle_p, NEWLINE);
651        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;  
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);
712                 /* tensor data, so have a 3x3 matrix to output as a row                  }
713                  * of 9 data points */              }
714                  int count = 0;              if (!Finley_noError())
715                  for (int m=0; m<shape; m++) {                  break;
716                    for (int n=0; n<shape; n++) {          }
717                       fprintf(fileHandle_p, " %e", sampleAvg[count]);      }
718                       count++;      /* only continue if no error occurred */
719                    }      if (writeCellData && Finley_noError()) {
720                    for (int n=0; n<3-shape; n++) fprintf(fileHandle_p, " %e", 0.);          strcat(txtBuffer, ">\n");
721                  }          if ( mpi_size > 1) {
722                  for (int m=0; m<3-shape; m++)              MPI_RANK0_WRITE_SHARED(txtBuffer);
723                     for (int n=0; n<3; n++) fprintf(fileHandle_p, " %e", 0.);          } else {
724                  }              fputs(txtBuffer, fileHandle_p);
725                fprintf(fileHandle_p, "\n");          }
726               }  
727               fprintf(fileHandle_p, "</DataArray>\n");          /* write the arrays */
728           }          for (dataIdx = 0; dataIdx < num_data; dataIdx++) {
729         }              if (!isEmpty(data_pp[dataIdx]) && isCellCentered[dataIdx]) {
730         fprintf(fileHandle_p, "</CellData>\n");                  dim_t numPointsPerSample=getNumDataPointsPerSample(data_pp[dataIdx]);
731    }                  dim_t rank = getDataPointRank(data_pp[dataIdx]);
732    /* point data */                  dim_t nComp = getDataPointSize(data_pp[dataIdx]);
733    if (write_pointdata) {                  dim_t nCompReqd = 1; /* number of components mpi_required */
734         /* mark the active data arrays */                  if (rank == 0) {
735         bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;                      nCompReqd = 1;
736         fprintf(fileHandle_p, "<PointData");                      shape = 0;
737         for (i_data =0 ;i_data<num_data;++i_data) {                  } else if (rank == 1) {
738              if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data]) {                      nCompReqd = 3;
739                  /* if the rank == 0:   --> scalar data                      shape = getDataPointShape(data_pp[dataIdx], 0);
740                   * if the rank == 1:   --> vector data                      if (shape > 3) {
741                   * if the rank == 2:   --> tensor data                          Finley_setError(VALUE_ERROR, "saveVTK: rank 1 objects must have 3 components at most.");
742                   */                      }
743                  switch(getDataPointRank(data_pp[i_data])) {                  } else {
744                     case 0:                      nCompReqd = 9;
745                         if (! set_scalar) {                      shape = getDataPointShape(data_pp[dataIdx], 0);
746                               fprintf(fileHandle_p," Scalars=\"%s\"",names_p[i_data]);                      if (shape > 3 || shape != getDataPointShape(data_pp[dataIdx], 1)) {
747                               set_scalar=TRUE;                          Finley_setError(VALUE_ERROR, "saveVTK: rank 2 objects of shape 2x2 or 3x3 supported only.");
748                         }                      }
749                         break;                  }
750                     case 1:                  /* bail out if an error occurred */
751                         if (! set_vector) {                  if (!Finley_noError())
752                               fprintf(fileHandle_p," Vectors=\"%s\"",names_p[i_data]);                      break;
753                               set_vector=TRUE;  
754                         }                  sprintf(txtBuffer, tag_Float_DataArray, names_p[dataIdx], nCompReqd);
755                         break;                  if ( mpi_size > 1) {
756                     case 2:                      MPI_RANK0_WRITE_SHARED(txtBuffer);
757                         if (! set_tensor) {                  } else {
758                               fprintf(fileHandle_p," Tensors=\"%s\"",names_p[i_data]);                      fputs(txtBuffer, fileHandle_p);
759                               set_tensor=TRUE;                  }
760                         }  
761                         break;                  txtBuffer[0] = '\0';
762                     default:                  txtBufferInUse = 0;
763                         sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);                  for (i=0; i<numCells; i++) {
764                         Finley_setError(VALUE_ERROR,error_msg);                      if (elements->Owner[i] == my_mpi_rank) {
765                         fclose(fileHandle_p);                          double *values = getSampleData(data_pp[dataIdx], i);
766                         return;                          for (l = 0; l < numCellFactor; l++) {
767                                double sampleAvg[NCOMP_MAX];
768                                dim_t nCompUsed = MIN(nComp, NCOMP_MAX);
769    
770                                /* average over number of points in the sample */
771                                if (isExpanded(data_pp[dataIdx])) {
772                                    dim_t hits=0, hits_old;
773                                    for (k=0; k<nCompUsed; k++) sampleAvg[k]=0;
774                                    for (j=0; j<numPointsPerSample; j++) {
775                                        hits_old=hits;
776                                        if (nodeInQuadrant(quadNodes_p, typeId, j, l)) {
777                                            hits++;
778                                            for (k=0; k<nCompUsed; k++) {
779                                                sampleAvg[k] += values[INDEX2(k,j,nComp)];
780                                            }
781                                        }
782                                    }
783                                    for (k=0; k<nCompUsed; k++)
784                                        sampleAvg[k] /= MAX(hits, 1);
785                                } else {
786                                    for (k=0; k<nCompUsed; k++)
787                                        sampleAvg[k] = values[k];
788                                } /* isExpanded */
789    
790                                /* if the number of required components is more than
791                                 * the number of actual components, pad with zeros
792                                 */
793                                /* probably only need to get shape of first element */
794                                if (nCompReqd == 1) {
795                                    sprintf(tmpBuffer, SCALAR_FORMAT, sampleAvg[0]);
796                                } else if (nCompReqd == 3) {
797                                    if (shape==1) {
798                                        sprintf(tmpBuffer, VECTOR_FORMAT,
799                                            sampleAvg[0], 0.f, 0.f);
800                                    } else if (shape==2) {
801                                        sprintf(tmpBuffer, VECTOR_FORMAT,
802                                            sampleAvg[0], sampleAvg[1], 0.f);
803                                    } else if (shape==3) {
804                                        sprintf(tmpBuffer, VECTOR_FORMAT,
805                                            sampleAvg[0],sampleAvg[1],sampleAvg[2]);
806                                    }
807                                } else if (nCompReqd == 9) {
808                                    if (shape==1) {
809                                        sprintf(tmpBuffer, TENSOR_FORMAT,
810                                            sampleAvg[0], 0.f, 0.f,
811                                                     0.f, 0.f, 0.f,
812                                                     0.f, 0.f, 0.f);
813                                    } else if (shape==2) {
814                                        sprintf(tmpBuffer, TENSOR_FORMAT,
815                                            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         fprintf(fileHandle_p, ">\n");                  break;
892         /* write the arrays */          }
893         for (i_data =0 ;i_data<num_data;++i_data) {      }
894            if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data]) {      /* only continue if no error occurred */
895               int numPointsPerSample = elements->ReferenceElement->numQuadNodes;      if (writePointData && Finley_noError()) {
896               int rank = getDataPointRank(data_pp[i_data]);          strcat(txtBuffer, ">\n");
897               int nComp = getDataPointSize(data_pp[i_data]);          if ( mpi_size > 1) {
898               int nCompReqd=1;   /* the number of components required by vtk */              MPI_RANK0_WRITE_SHARED(txtBuffer);
899               int shape=0;          } else {
900               if (rank == 0) {              fputs(txtBuffer, fileHandle_p);
901                  nCompReqd = 1;          }
902               } else if (rank == 1) {  
903                   shape=getDataPointShape(data_pp[i_data], 0);          /* write the arrays */
904                   if  (shape>3) {          for (dataIdx = 0; dataIdx < num_data; dataIdx++) {
905                       Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");              if (!isEmpty(data_pp[dataIdx]) && !isCellCentered[dataIdx]) {
906                       fclose(fileHandle_p);                  Finley_NodeMapping* nodeMapping;
907                       return;                  dim_t rank = getDataPointRank(data_pp[dataIdx]);
908                   }                  dim_t nCompReqd = 1; /* number of components mpi_required */
909                   nCompReqd = 3;                  if (getFunctionSpaceType(data_pp[dataIdx]) == FINLEY_REDUCED_NODES) {
910               } else {                      nodeMapping = mesh_p->Nodes->reducedNodesMapping;
911                   shape=getDataPointShape(data_pp[i_data], 0);                  } else {
912                   if  (shape>3 || shape != getDataPointShape(data_pp[i_data], 1)) {                      nodeMapping = mesh_p->Nodes->nodesMapping;
913                       Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");                  }
914                       fclose(fileHandle_p);                  if (rank == 0) {
915                       return;                      nCompReqd = 1;
916                   }                      shape = 0;
917                   nCompReqd = 9;                  } else if (rank == 1) {
918               }                      nCompReqd = 3;
919               fprintf(fileHandle_p, "<DataArray Name=\"%s\" type=\"Float32\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);                      shape = getDataPointShape(data_pp[dataIdx], 0);
920           /* write out the data */                      if (shape > 3) {
921           /* if the number of required components is more than the number                          Finley_setError(VALUE_ERROR, "saveVTK: rank 1 objects must have 3 components at most.");
922            * of actual components, pad with zeros                      }
923                 */                  } else {
924               bool_t do_write=TRUE;                      nCompReqd = 9;
925               for (i=0; i<mesh_p->Nodes->numNodes; i++) {                      shape=getDataPointShape(data_pp[dataIdx], 0);
926                 if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {                      if (shape > 3 || shape != getDataPointShape(data_pp[dataIdx], 1)) {
927                  if (mesh_p->Nodes->toReduced[i]>=0) {                          Finley_setError(VALUE_ERROR, "saveVTK: rank 2 objects of shape 2x2 or 3x3 supported only.");
928                         switch(getFunctionSpaceType(data_pp[i_data])) {                      }
929                            case FINLEY_DEGREES_OF_FREEDOM:                  }
930                               values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);                  /* bail out if an error occurred */
931                               break;                  if (!Finley_noError())
932                            case FINLEY_REDUCED_DEGREES_OF_FREEDOM:                      break;
933                               values = getSampleData(data_pp[i_data],mesh_p->Nodes->reducedDegreeOfFreedom[i]);  
934                               break;                  sprintf(txtBuffer, tag_Float_DataArray, names_p[dataIdx], nCompReqd);
935                            case FINLEY_NODES:                  if ( mpi_size > 1) {
936                           values = getSampleData(data_pp[i_data],i);                      MPI_RANK0_WRITE_SHARED(txtBuffer);
937                               break;                  } else {
938                         }                      fputs(txtBuffer, fileHandle_p);
939                         do_write=TRUE;                  }
940                      } else {  
941                         do_write=FALSE;                  txtBuffer[0] = '\0';
942                      }                  txtBufferInUse = 0;
943                 } else {                  for (i=0; i<mesh_p->Nodes->numNodes; i++) {
944                      do_write=TRUE;                      k = globalNodeIndex[i];
945                      switch(getFunctionSpaceType(data_pp[i_data])) {                      if ( (myFirstNode <= k) && (k < myLastNode) ) {
946                         case FINLEY_DEGREES_OF_FREEDOM:                          double *values = getSampleData(data_pp[dataIdx], nodeMapping->target[i]);
947                        values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);                          /* if the number of mpi_required components is more than
948                            break;                           * the number of actual components, pad with zeros.
949                         case FINLEY_NODES:                           * Probably only need to get shape of first element */
950                        values = getSampleData(data_pp[i_data],i);                          if (nCompReqd == 1) {
951                            break;                              sprintf(tmpBuffer, SCALAR_FORMAT, values[0]);
952                      }                          } else if (nCompReqd == 3) {
953                 }                              if (shape==1) {
954             /* write the data different ways for scalar, vector and tensor */                                  sprintf(tmpBuffer, VECTOR_FORMAT,
955                 if (do_write) {                                          values[0], 0.f, 0.f);
956                if (nCompReqd == 1) {                              } else if (shape==2) {
957                  fprintf(fileHandle_p, " %e", values[0]);                                  sprintf(tmpBuffer, VECTOR_FORMAT,
958                } else if (nCompReqd == 3) {                                          values[0], values[1], 0.f);
959                  for (int m=0; m<shape; m++) fprintf(fileHandle_p, " %e", values[m]);                              } else if (shape==3) {
960                  for (int m=0; m<nCompReqd-shape; m++) fprintf(fileHandle_p, " %e", 0.);                                  sprintf(tmpBuffer, VECTOR_FORMAT,
961                } else if (nCompReqd == 9) {                                          values[0], values[1], values[2]);
962                  /* tensor data, so have a 3x3 matrix to output as a row                              }
963                   * of 9 data points */                          } else if (nCompReqd == 9) {
964                  int count = 0;                              if (shape==1) {
965                  for (int m=0; m<shape; m++) {                                  sprintf(tmpBuffer, TENSOR_FORMAT,
966                    for (int n=0; n<shape; n++) {                                      values[0], 0.f, 0.f,
967                      fprintf(fileHandle_p, " %e", values[count]);                                            0.f, 0.f, 0.f,
968                      count++;                                            0.f, 0.f, 0.f);
969                    }                              } else if (shape==2) {
970                    for (int n=0; n<3-shape; n++) fprintf(fileHandle_p, " %e", 0.);                                  sprintf(tmpBuffer, TENSOR_FORMAT,
971                  }                                      values[0], values[1], 0.f,
972                  for (int m=0; m<3-shape; m++)                                        values[2], values[3], 0.f,
973                      for (int n=0; n<3; n++) fprintf(fileHandle_p, " %e", 0.);                                            0.f,       0.f, 0.f);
974                }                              } else if (shape==3) {
975                    fprintf(fileHandle_p, "\n");                                  sprintf(tmpBuffer, TENSOR_FORMAT,
976                 }                                      values[0], values[1], values[2],
977               }                                      values[3], values[4], values[5],
978               fprintf(fileHandle_p, "</DataArray>\n");                                      values[6], values[7], values[8]);
979             }                              }
980         }                          }
981         fprintf(fileHandle_p, "</PointData>\n");                          if ( mpi_size > 1) {
982    }                              __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
983    /* finish off the piece */                          } else {
984    fprintf(fileHandle_p, "</Piece>\n");                              fputs(tmpBuffer, fileHandle_p);
985                            }
986    fprintf(fileHandle_p, "</UnstructuredGrid>\n");                      } /* if this is my node */
987    /* write the xml footer */                  } /* for i (numNodes) */
988    fprintf(fileHandle_p, "</VTKFile>\n");  
989    /* close the file */                  if ( mpi_size > 1) {
990    fclose(fileHandle_p);                      MPI_WRITE_ORDERED(txtBuffer, txtBufferInUse);
991    return;                      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.201  
changed lines
  Added in v.2141

  ViewVC Help
Powered by ViewVC 1.1.26