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

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

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

trunk/esys2/finley/src/finleyC/Mesh_saveVTK.c revision 121 by jgs, Fri May 6 04:26:16 2005 UTC trunk/finley/src/Mesh_saveVTK.c revision 796 by dhawcroft, Mon Jul 31 06:16:17 2006 UTC
# Line 1  Line 1 
1  /* $Id$ */  /*
2    ************************************************************
3    *          Copyright 2006 by ACcESS MNRF                   *
4    *                                                          *
5    *              http://www.access.edu.au                    *
6    *       Primary Business: Queensland, Australia            *
7    *  Licensed under the Open Software License version 3.0    *
8    *     http://www.opensource.org/licenses/osl-3.0.php       *
9    *                                                          *
10    ************************************************************
11    */
12    
13    
14  /**************************************************************/  /**************************************************************/
15    
# Line 6  Line 17 
17    
18  /**************************************************************/  /**************************************************************/
19    
 /*   Copyrights by ACcESS, Australia 2004 */  
20  /*   Author: Paul Cochrane, cochrane@esscc.uq.edu.au */  /*   Author: Paul Cochrane, cochrane@esscc.uq.edu.au */
21    /*   MPI-IO version: Derick Hawcroft, d.hawcroft@uq.edu.au         */
22    
23    /*   Version: $Id$ */
24    
25  /**************************************************************/  /**************************************************************/
26    
27  #include "Finley.h"  
 #include "Common.h"  
28  #include "Mesh.h"  #include "Mesh.h"
 #include "escript/Data/DataC.h"  
29  #include "vtkCellType.h"  /* copied from vtk source directory !!! */  #include "vtkCellType.h"  /* copied from vtk source directory !!! */
30    
31  void Finley_Mesh_saveVTK(const char * filename_p, Finley_Mesh *mesh_p, escriptDataC* data_p) {  /*
32    /* if there is no mesh we just return */   MPI version notes:
33    if (mesh_p==NULL) return;  
34    Finley_ElementFile* elements=NULL;   ******************************************************************************
35    char elemTypeStr[32];   ***                                                                       ****
36    int i, j, k, numVTKNodesPerElement, isCellCentered, nodetype;   *** WARNING: Won't work for meshes with peridodic boundary conditions yet ****
37     ***                                                                       ****  
38     ******************************************************************************
39    
40     In this version, the rank==0 process writes *all* opening and closing
41     XML tags.
42     Individual process data is copied to a buffer before being written
43     out. The  routines are collectively called and will be called in the natural
44     ordering i.e 0 to maxProcs-1.
45    
46    */
47    
48    #ifdef PASO_MPI
49    
50    
51    //#define MPIO_HINTS
52    
53    
54    
55    #define MPIO_DEBUG(str) \
56    { \
57        if(myRank == 0) \
58        printf("==== MPI-IO => %s \n", str); \
59    }
60    
61    void Finley_Mesh_saveVTK_MPIO(const char * filename_p, Finley_Mesh *mesh_p, const dim_t num_data,char* *names_p, escriptDataC* *data_pp)
62    {
63      double time0 = Paso_timer();
64      int    numPoints,
65      numCells = -1,
66                 myRank,comm,gsize,
67                 numLocal,
68                 nDim,
69                 shape;
70      size_t __n;
71      int* failSend;
72      int i,j,k,m,n,count;
73      int numGlobalCells = 0;
74      index_t  *nodesGlobal=NULL;   // used to get the connectivity  right for VTK
75    
76      /* variables associatted with write_celldata/pointdata */
77      int numPointsPerSample,
78      nComp,
79      nCompReqd;
80    double* values, rtmp;    double* values, rtmp;
   int nDim = mesh_p->Nodes->numDim;  
81    
82    /* get a pointer to the relevant mesh component */    // Local element info (for debugging)
83    if (isEmpty(data_p)) {    size_t numLocalCells,
84      elements=mesh_p->Elements;    numInternalCells,
85    } else {    numBoundaryCells;
86      switch(getFunctionSpaceType(data_p)) {  
87      case(FINLEY_DEGREES_OF_FREEDOM):    int rank;
88        nodetype = FINLEY_DEGREES_OF_FREEDOM;  
89        isCellCentered = FALSE;    int amode = MPI_MODE_CREATE | MPI_MODE_WRONLY |  MPI_MODE_SEQUENTIAL;
90        elements = mesh_p->Elements;  
91        break;    comm   = mesh_p->Nodes->MPIInfo->comm;
92      case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):    myRank = mesh_p->Nodes->MPIInfo->rank;
93        Finley_ErrorCode=VALUE_ERROR;    gsize  = mesh_p->Nodes->MPIInfo->size;
94        sprintf(Finley_ErrorMsg,  
95            "Reduced degrees of freedom is not yet "    MPI_File fh;
96            "implemented for saving vtk files\n");    MPI_Status status;
97        return;    MPI_Request req;
98      case(FINLEY_NODES):    MPI_Info infoHints;
99        nodetype=FINLEY_NODES;  
100        isCellCentered=FALSE;    char error_msg[LenErrorMsg_MAX];
101        elements=mesh_p->Elements;  
102        break;    int i_data;
103      case(FINLEY_ELEMENTS):  
104        isCellCentered=TRUE;    int nodetype=FINLEY_DEGREES_OF_FREEDOM;
105        elements=mesh_p->Elements;    int elementtype=FINLEY_UNKNOWN;
106        break;    bool_t isCellCentered[num_data],write_celldata=FALSE,write_pointdata=FALSE;
107      case(FINLEY_FACE_ELEMENTS):  
108        isCellCentered=TRUE;    ElementTypeId TypeId;
109        elements=mesh_p->FaceElements;  
110        break;    int numVTKNodesPerElement;
111      case(FINLEY_POINTS):    int cellType;
112        isCellCentered=TRUE;    char elemTypeStr[32];
113        elements=mesh_p->Points;  
114        break;    Finley_ElementFile* elements=NULL;
115      case(FINLEY_CONTACT_ELEMENTS_1):  
116      case(FINLEY_CONTACT_ELEMENTS_2):  
117        isCellCentered=TRUE;    // Local node info
118        elements=mesh_p->ContactElements;    int numInternalNodes,
119        break;    numLocalNodes,
120      default:    numBoundaryNodes,
121        Finley_ErrorCode=TYPE_ERROR;    localDOF;
122        sprintf(Finley_ErrorMsg,  
123            "Finley does not know anything about function space type %d",  
124            getFunctionSpaceType(data_p));    nDim  = mesh_p->Nodes->numDim;
125        return;  
126    #ifdef MPIO_HINTS
127      MPI_Info_create(&infoHints);
128      //  MPI_Info_set(infoHints, "striping_unit",        "424288");
129      //  MPI_Info_set(infoHints, "striping_factor",      "16");
130      //  MPI_Info_set(infoHints, "collective_buffering", "true");
131      //  MPI_Info_set(infoHints, "cb_block_size",        "131072");
132      //  MPI_Info_set(infoHints, "cb_buffer_size",       "1048567");
133      //  MPI_Info_set(infoHints, "cb_nodes",             "8");
134      //    MPI_Info_set(infoHints, "access_style", "write_once, sequential");
135    
136      //XFS only
137      //   MPI_Info_set(infoHints, "direct_write",          "true");
138    #else
139      infoHints = MPI_INFO_NULL;
140    #endif
141    
142      // Holds a local node/element index to their global arrays
143      struct localIndexCache
144      {
145        index_t *values;
146        int size;
147      };
148      typedef struct localIndexCache localIndexCache;
149    
150      localIndexCache nodeCache,
151      elementCache;
152    
153      // Collective Call
154      MPI_File_open(mesh_p->Nodes->MPIInfo->comm, (char*)filename_p, amode,infoHints, &fh);
155      MPI_File_set_view(fh,MPI_DISPLACEMENT_CURRENT,MPI_CHAR, MPI_CHAR, "native" , infoHints);
156    
157      MPIO_DEBUG(" ***** Enter saveVTK ******")
158    
159      for (i_data=0;i_data<num_data;++i_data)
160      {
161        if (! isEmpty(data_pp[i_data]) )
162        {
163          switch(getFunctionSpaceType(data_pp[i_data]))
164          {
165          case FINLEY_DEGREES_OF_FREEDOM:
166            nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
167            if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)
168            {
169              elementtype=FINLEY_ELEMENTS;
170            }
171            else
172            {
173              Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
174              return ;
175            }
176            isCellCentered[i_data]=FALSE;
177            break;
178          case FINLEY_REDUCED_DEGREES_OF_FREEDOM:
179            nodetype = FINLEY_REDUCED_DEGREES_OF_FREEDOM;
180            if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)
181            {
182              elementtype=FINLEY_ELEMENTS;
183            }
184            else
185            {
186              Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
187              return;
188            }
189            isCellCentered[i_data]=FALSE;
190            break;
191          case FINLEY_NODES:
192            nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
193            if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)
194            {
195              elementtype=FINLEY_ELEMENTS;
196            }
197            else
198            {
199              Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
200              return;
201            }
202            isCellCentered[i_data]=FALSE;
203            break;
204          case FINLEY_ELEMENTS:
205            nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
206            if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)
207            {
208              elementtype=FINLEY_ELEMENTS;
209            }
210            else
211            {
212              Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
213              return;
214            }
215            isCellCentered[i_data]=TRUE;
216            break;
217          case FINLEY_FACE_ELEMENTS:
218            nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
219            if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_FACE_ELEMENTS)
220            {
221              elementtype=FINLEY_FACE_ELEMENTS;
222            }
223            else
224            {
225              Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
226              return;
227    
228            }
229            isCellCentered[i_data]=TRUE;
230            break;
231          case FINLEY_POINTS:
232            nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
233            if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_POINTS)
234            {
235              elementtype=FINLEY_POINTS;
236            }
237            else
238            {
239              Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
240              return;
241    
242            }
243            isCellCentered[i_data]=TRUE;
244            break;
245          case FINLEY_CONTACT_ELEMENTS_1:
246            nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
247            if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1)
248            {
249              elementtype=FINLEY_CONTACT_ELEMENTS_1;
250            }
251            else
252            {
253              Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
254              return;
255    
256            }
257            isCellCentered[i_data]=TRUE;
258            break;
259          case FINLEY_CONTACT_ELEMENTS_2:
260            nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
261            if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1)
262            {
263              elementtype=FINLEY_CONTACT_ELEMENTS_1;
264            }
265            else
266            {
267              Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
268              return;
269    
270            }
271            isCellCentered[i_data]=TRUE;
272            break;
273          default:
274            sprintf(error_msg,"saveVTK: Finley does not know anything about function space type %d",getFunctionSpaceType(data_pp[i_data]));
275            Finley_setError(TYPE_ERROR,error_msg);
276            return;
277    
278          }
279    
280          if (isCellCentered[i_data])
281          {
282            write_celldata=TRUE;
283          }
284          else
285          {
286            write_pointdata=TRUE;
287          }
288      }      }
289    }    }
290    
291    /* the number of points */    Finley_NodeDistribution *dist;
292    int numPoints = mesh_p->Nodes->numNodes;    if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
293      {
294    /* the number of cells */      dist = mesh_p->Nodes->reducedDegreeOfFreedomDistribution;
295    if (elements == NULL) {    }
296      Finley_ErrorCode = VALUE_ERROR;    else
297      sprintf(Finley_ErrorMsg,    {
298          "elements object is NULL; cannot proceed");      dist = mesh_p->Nodes->degreeOfFreedomDistribution;
299      }
300    
301      numInternalNodes = dist->numInternal;
302      numBoundaryNodes = dist->numBoundary;
303    
304      localDOF =  dist->numLocal;
305    
306      numPoints        = dist->numGlobal;
307    
308      if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;
309      switch(elementtype)
310      {
311      case FINLEY_ELEMENTS:
312        elements=mesh_p->Elements;
313        break;
314      case FINLEY_FACE_ELEMENTS:
315        elements=mesh_p->FaceElements;
316        break;
317      case FINLEY_POINTS:
318        elements=mesh_p->Points;
319        break;
320      case FINLEY_CONTACT_ELEMENTS_1:
321        elements=mesh_p->ContactElements;
322        break;
323      }
324      if (elements==NULL)
325      {
326        Finley_setError(SYSTEM_ERROR,"saveVTK: undefined element file");
327        return ;
328      }
329    
330      numCells =  elements->numElements;
331      numGlobalCells = elements->elementDistribution->vtxdist[gsize];
332      numLocalCells    = elements->elementDistribution->numLocal;
333      numInternalCells = elements->elementDistribution->numInternal;
334      numBoundaryCells = elements->elementDistribution->numBoundary;
335    
336      if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
337      {
338        TypeId = elements->LinearReferenceElement->Type->TypeId;
339      }
340      else
341      {
342        TypeId = elements->ReferenceElement->Type->TypeId;
343      }
344    
345      switch(TypeId)
346      {
347      case Point1:
348      case Line2Face:
349      case Line3Face:
350      case Point1_Contact:
351      case Line2Face_Contact:
352      case Line3Face_Contact:
353        cellType = VTK_VERTEX;
354        numVTKNodesPerElement = 1;
355        strcpy(elemTypeStr, "VTK_VERTEX");
356        break;
357    
358      case Line2:
359      case Tri3Face:
360      case Rec4Face:
361      case Line2_Contact:
362      case Tri3_Contact:
363      case Tri3Face_Contact:
364      case Rec4Face_Contact:
365        cellType = VTK_LINE;
366        numVTKNodesPerElement = 2;
367        strcpy(elemTypeStr, "VTK_LINE");
368        break;
369    
370      case Tri3:
371      case Tet4Face:
372      case Tet4Face_Contact:
373        cellType = VTK_TRIANGLE;
374        numVTKNodesPerElement = 3;
375        strcpy(elemTypeStr, "VTK_TRIANGLE");
376        break;
377    
378      case Rec4:
379      case Hex8Face:
380      case Rec4_Contact:
381      case Hex8Face_Contact:
382        cellType = VTK_QUAD;
383        numVTKNodesPerElement = 4;
384        strcpy(elemTypeStr, "VTK_QUAD");
385        break;
386    
387      case Tet4:
388        cellType = VTK_TETRA;
389        numVTKNodesPerElement = 4;
390        strcpy(elemTypeStr, "VTK_TETRA");
391        break;
392    
393      case Hex8:
394        cellType = VTK_HEXAHEDRON;
395        numVTKNodesPerElement = 8;
396        strcpy(elemTypeStr, "VTK_HEXAHEDRON");
397        break;
398    
399      case Line3:
400      case Tri6Face:
401      case Rec8Face:
402      case Line3_Contact:
403      case Tri6Face_Contact:
404      case Rec8Face_Contact:
405        cellType = VTK_QUADRATIC_EDGE;
406        numVTKNodesPerElement = 3;
407        strcpy(elemTypeStr, "VTK_QUADRATIC_EDGE");
408        break;
409    
410      case Tri6:
411      case Tet10Face:
412      case Tri6_Contact:
413      case Tet10Face_Contact:
414        cellType = VTK_QUADRATIC_TRIANGLE;
415        numVTKNodesPerElement = 6;
416        strcpy(elemTypeStr, "VTK_QUADRATIC_TRIANGLE");
417        break;
418    
419      case Rec8:
420      case Hex20Face:
421      case Rec8_Contact:
422      case Hex20Face_Contact:
423        cellType = VTK_QUADRATIC_QUAD;
424        numVTKNodesPerElement = 8;
425        strcpy(elemTypeStr, "VTK_QUADRATIC_QUAD");
426        break;
427    
428      case Tet10:
429        cellType = VTK_QUADRATIC_TETRA;
430        numVTKNodesPerElement = 10;
431        strcpy(elemTypeStr, "VTK_QUADRATIC_TETRA");
432        break;
433    
434      case Hex20:
435        cellType = VTK_QUADRATIC_HEXAHEDRON;
436        numVTKNodesPerElement = 20;
437        strcpy(elemTypeStr, "VTK_QUADRATIC_HEXAHEDRON");
438        break;
439    
440      default:
441        sprintf(error_msg, "saveVTK: Element type %s is not supported by VTK",elements->ReferenceElement->Type->Name);
442        Finley_setError(VALUE_ERROR,error_msg);
443      return;      return;
444    }    }
445    int numCells = elements->numElements;    
446        /* Write XML Header */
447      if(myRank == 0)
448      {
449        char header[400];
450    
451        sprintf(header,"<?xml version=\"1.0\"?>\n" \
452                "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\">\n" \
453                "<UnstructuredGrid>\n" \
454                "<Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n" \
455                "<Points>\n" \
456                "<DataArray NumberOfComponents=\"%d\" type=\"Float32\" format=\"ascii\">\n"
457                ,numPoints,numGlobalCells,MAX(3,nDim));
458    
459    
460        MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);
461        MPI_Wait(&req,&status);
462      }
463    
464      MPIO_DEBUG(" Writing Coordinate Points... ")
465    
466      numLocalNodes=localDOF;
467    
468      //  values vary from 13-14 chars hence the strlen()
469      char* largebuf = MEMALLOC( numLocalNodes*14*nDim + numLocalNodes*2 + 1 ,char);
470      largebuf[0] = '\0';
471      char tmpbuf[14];
472      int tsz=0;
473      int numNodesOutput=0;
474      index_t pos=0;
475    
476      index_t *vtxdist = NULL, *DOFNodes=NULL,*forwardBuffer=NULL,*backwardBuffer=NULL;
477    
478      DOFNodes   = MEMALLOC(mesh_p->Nodes->numNodes,index_t);
479      nodeCache.values = MEMALLOC( numLocalNodes, index_t);
480      index_t bc_pos = 0;
481    
482      // Custom strcat:  avoids expensive strlen(3) call by  strcat(3)
483      int __len,__j;
484      char  *zero = "0.000000e+00 ";
485      char  *newline = "\n";
486    #define __STRCAT(__buf,x)  \
487    {                  \
488       __j = -1;      \
489       while(__j++ < __len)  \
490        *(__buf+tsz+__j)=*(x+__j); \
491    }
492    
493      for (i = 0; i < mesh_p->Nodes->numNodes; i++)
494      {
495        // This is the bit that will break for periodic BCs because it assumes that there is a one to one
496        // correspondance between nodes and Degrees of freedom
497        DOFNodes[mesh_p->Nodes->degreeOfFreedom[i]] = i;
498    
499        /* local node ?*/
500        if( mesh_p->Nodes->degreeOfFreedom[i] < localDOF )
501        {
502          for (j = 0; j < nDim; j++)
503          {
504            sprintf(tmpbuf,"%e ", mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)] );
505            __len = strlen(tmpbuf);
506            __STRCAT(largebuf,tmpbuf)
507            tsz+=__len;
508          }
509          for (k=0; k<3-nDim; k++)
510          {
511            __len = 13;
512            __STRCAT(largebuf,zero)
513            tsz+=__len;
514          }
515          __len=1;
516          __STRCAT(largebuf,newline)
517          tsz += 1;
518          nodeCache.values[numNodesOutput++]=i;
519        }
520      }
521    
522      nodeCache.size=numNodesOutput;
523    
524      largebuf[tsz] = '\0';
525      MPI_File_write_ordered(fh, largebuf,tsz, MPI_CHAR, &status);
526      MEMFREE(largebuf);
527    
528      nodesGlobal = MEMALLOC(mesh_p->Nodes->numNodes ,index_t);
529    
530      // form distribution info on who output which nodes
531      vtxdist = MEMALLOC( gsize+1, index_t );
532      vtxdist[0]=0;
533      MPI_Allgather(&numNodesOutput,1,MPI_INT,vtxdist+1,1,MPI_INT,comm);
534      for( i=0; i<gsize; i++ )
535        vtxdist[i+1]+=vtxdist[i];
536    
537      // will not work for periodic boundary conditions
538      // calculate the local nodes file positions
539      pos = 0;
540      for( i=0; i<mesh_p->Nodes->numNodes; i++ )
541      {
542        if( mesh_p->Nodes->degreeOfFreedom[i]< localDOF )
543        {
544          nodesGlobal[i] = vtxdist[myRank] + pos++;
545        }
546        else
547          nodesGlobal[i] = -1;
548      }
549    
550      // communicate the local Nodes file position to the interested parties
551      // send local info
552      forwardBuffer = MEMALLOC( mesh_p->Nodes->numNodes, index_t );
553      for( n=0; n < dist->numNeighbours; n++ )
554      {
555        if(  dist->edges[n]->numForward)
556        {
557          for( i=0; i < dist->edges[n]->numForward; i++ )
558            forwardBuffer[i] = nodesGlobal[DOFNodes[dist->edges[n]->indexForward[i] ]];
559          Paso_CommBuffer_pack( mesh_p->Nodes->CommBuffer, dist->neighbours[n], NULL, forwardBuffer, sizeof(index_t), 0 );
560          Paso_CommBuffer_send( mesh_p->Nodes->CommBuffer, dist->neighbours[n], sizeof(index_t) );
561        }
562      }
563      // receive external info
564      backwardBuffer = MEMALLOC( mesh_p->Nodes->numNodes, index_t );
565      for( n=0; n < dist->numNeighbours; n++ )
566      {
567        if( dist->edges[n]->numBackward )
568        {
569          Paso_CommBuffer_recv(mesh_p->Nodes->CommBuffer, dist->neighbours[n], sizeof(index_t));
570          Paso_CommBuffer_unpack(mesh_p->Nodes->CommBuffer, dist->neighbours[n], NULL, backwardBuffer, sizeof(index_t), 0 );
571          for( i=0; i<dist->edges[n]->numBackward; i++ )
572            nodesGlobal[DOFNodes[dist->edges[n]->indexBackward[i] ]] = backwardBuffer[i];
573        }
574      }
575    
576      MEMFREE(vtxdist);
577      MEMFREE(DOFNodes);
578      MEMFREE(backwardBuffer);
579      MEMFREE(forwardBuffer);
580    
581      if( myRank == 0)
582      {
583        char* tags = "</DataArray>\n</Points>\n<Cells>\n<DataArray Name=\"connectivity\" type=\"Int32\" " \
584                     "format=\"ascii\">\n" ;
585        MPI_File_iwrite_shared(fh,tags,strlen(tags),MPI_CHAR,&req);
586        MPI_Wait(&req,&status);
587      }
588      MPIO_DEBUG(" Done Writing Coordinate Points ")
589    
590      /* BEGIN CONNECTIVITY */
591    
592      int NN = elements->ReferenceElement->Type->numNodes; /* num Nodes holding ref-element */
593    
594      // Collective
595      MPIO_DEBUG(" Writing Connectivity... ")
596    
597      // TODO: Improve on upper bound , will fail for very large meshes!!
598      size_t sz = numLocalCells*6*numVTKNodesPerElement + numLocalCells;
599      largebuf = MEMALLOC(sz,char);
600      largebuf[0] = '\0';
601      tsz=0;
602      pos = 0;
603      // numCells?
604      elementCache.values = MEMALLOC(numLocalCells,index_t);
605      if (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM)
606      {
607        for (i = 0; i < numCells; i++)
608        {
609    
610          if (elements->Id[i] >= elements->elementDistribution->vtxdist[i] &&  elements->Id[i] <= elements->elementDistribution->vtxdist[i+1] - 1 )
611          {
612            for (j = 0; j < numVTKNodesPerElement; j++)
613            {
614              sprintf(tmpbuf,"%d ",nodesGlobal[mesh_p->Nodes->toReduced[elements->Nodes[INDEX2(elements->ReferenceElement->Type->linearNodes[j], i, NN)]]]);
615              __len=strlen(tmpbuf);
616              __STRCAT(largebuf,tmpbuf)
617              tsz+=__len;
618            }
619            __len=1;
620            __STRCAT(largebuf,newline)
621            tsz+=1;
622    
623            elementCache.values[pos++]=i;
624          }
625        }
626      }
627      else if (VTK_QUADRATIC_HEXAHEDRON==cellType)
628      {
629        char tmpbuf2[20*20];
630    
631        for (i = 0; i < numCells; i++)
632        {
633          if( elements->Id[i] >= elements->elementDistribution->vtxdist[myRank] && elements->Id[i] <= elements->elementDistribution->vtxdist[myRank+1]-1)
634          {
635            sprintf(tmpbuf2,"%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d\n",
636                    nodesGlobal[elements->Nodes[INDEX2(0, i, NN)]],
637                    nodesGlobal[elements->Nodes[INDEX2(1, i, NN)]],
638                    nodesGlobal[elements->Nodes[INDEX2(2, i, NN)]],
639                    nodesGlobal[elements->Nodes[INDEX2(3, i, NN)]],
640                    nodesGlobal[elements->Nodes[INDEX2(4, i, NN)]],
641                    nodesGlobal[elements->Nodes[INDEX2(5, i, NN)]],
642                    nodesGlobal[elements->Nodes[INDEX2(6, i, NN)]],
643                    nodesGlobal[elements->Nodes[INDEX2(7, i, NN)]],
644                    nodesGlobal[elements->Nodes[INDEX2(8, i, NN)]],
645                    nodesGlobal[elements->Nodes[INDEX2(9, i, NN)]],
646                    nodesGlobal[elements->Nodes[INDEX2(10, i, NN)]],
647                    nodesGlobal[elements->Nodes[INDEX2(11, i, NN)]],
648                    nodesGlobal[elements->Nodes[INDEX2(16, i, NN)]],
649                    nodesGlobal[elements->Nodes[INDEX2(17, i, NN)]],
650                    nodesGlobal[elements->Nodes[INDEX2(18, i, NN)]],
651                    nodesGlobal[elements->Nodes[INDEX2(19, i, NN)]],
652                    nodesGlobal[elements->Nodes[INDEX2(12, i, NN)]],
653                    nodesGlobal[elements->Nodes[INDEX2(13, i, NN)]],
654                    nodesGlobal[elements->Nodes[INDEX2(14, i, NN)]],
655                    nodesGlobal[elements->Nodes[INDEX2(15, i, NN)]]);
656            __len=strlen(tmpbuf2);
657            __STRCAT(largebuf,tmpbuf2)
658            tsz+=__len;
659            elementCache.values[pos++]=i;
660          }
661        }
662      }
663      else if (numVTKNodesPerElement!=NN)
664      {
665    
666        for (i = 0; i < numCells; i++)
667        {
668          for (j = 0; j < numVTKNodesPerElement; j++)
669          {
670            sprintf(tmpbuf,"%d ", nodesGlobal[elements->Nodes[INDEX2(elements->ReferenceElement->Type->geoNodes[j], i, NN)]]);
671            __len=strlen(tmpbuf);
672            __STRCAT(largebuf,tmpbuf)
673            tsz+=__len;
674          }
675          __len=1;
676          __STRCAT(largebuf,newline)
677          tsz+=1;
678          elementCache.values[pos++]=i;
679        }
680      }
681      else
682      {
683    
684        for(i = 0;i  < numCells ; i++)
685        {
686          // is this element in domain of process with "myRank"
687    
688          if( elements->Id[i] >= elements->elementDistribution->vtxdist[myRank] && elements->Id[i] <= elements->elementDistribution->vtxdist[myRank+1]-1)
689          {
690            for (j = 0; j < numVTKNodesPerElement; j++)
691            {
692              sprintf(tmpbuf,"%d ", nodesGlobal[ elements->Nodes[INDEX2(j, i, NN) ] ] );
693              __len=strlen(tmpbuf);
694              __STRCAT(largebuf,tmpbuf)
695              tsz+=__len;
696            }
697            __len=1;
698            __STRCAT(largebuf,newline)
699            tsz+=1;
700            elementCache.values[pos++]=i;
701          }
702        }
703      }
704      elementCache.size = pos;
705    
706      MPI_File_write_ordered(fh,largebuf,tsz, MPI_CHAR, &status);
707      MEMFREE(largebuf);
708      MPIO_DEBUG(" Done Writing Connectivity ")
709      MPIO_DEBUG(" Writing Offsets & Types... ")
710    
711      // Non-Collective
712      if( myRank == 0)
713      {
714        // write out the DataArray element for the offsets
715        char* tag1 = "</DataArray>\n<DataArray Name=\"offsets\" type=\"Int32\" format=\"ascii\">\n";
716        char* tag2 = "</DataArray>\n";
717        char *tag3 =  "<DataArray Name=\"types\" type=\"UInt8\" format=\"ascii\">\n";
718        char *tag4 = "</DataArray>\n</Cells>\n";
719    
720        int n = numVTKNodesPerElement;
721    
722        int sz=0;
723        int lg = log10(numGlobalCells * n) + 1;
724        sz += numGlobalCells*lg;
725        sz += numGlobalCells;
726        tsz = 0;
727    
728        char* largebuf = MEMALLOC(sz  + strlen(tag1) + strlen(tag2) + strlen(tag3) + strlen(tag4),char);
729        largebuf[0] ='\0';
730        char tmp[10];
731    
732        __len = strlen(tag1);
733        __STRCAT(largebuf,tag1)
734        tsz +=  __len;
735    
736        for (i=numVTKNodesPerElement; i<=numGlobalCells*numVTKNodesPerElement; i+=numVTKNodesPerElement)
737        {
738          sprintf(tmp,"%d\n", i);
739          __len=strlen(tmp);
740          __STRCAT(largebuf,tmp)
741          tsz+=__len;
742        }
743    
744        __len=strlen(tag2);
745        __STRCAT(largebuf,tag2)
746        tsz+=__len;
747    
748        MPI_File_iwrite_shared(fh,largebuf, tsz,MPI_CHAR,&req);
749        MPI_Wait(&req,&status);
750        
751        // Re-using buffer!!
752        largebuf[0] = '\0';
753        tsz = 0;
754        __len = strlen(tag3);
755        __STRCAT(largebuf,tag3)
756        tsz+=__len;
757        for (i=0; i<numGlobalCells; i++)
758        {
759          sprintf(tmp, "%d\n", cellType);
760          __len=strlen(tmp);
761          __STRCAT(largebuf,tmp)
762          tsz+=__len;
763        }
764        __len=strlen(tag4);
765        __STRCAT(largebuf,tag4)
766        tsz+=__len;
767    
768        MPI_File_iwrite_shared(fh,largebuf,tsz,MPI_CHAR,&req);
769        MPI_Wait(&req,&status);
770        MEMFREE(largebuf);
771      }
772    
773      MPIO_DEBUG(" Done Writing Offsets & Types ")
774    
775      // Write Point Data Header Tags
776      if( myRank == 0)
777      {
778        char header[600];
779        char tmpBuf[50];
780    
781        if (write_pointdata)
782        {
783          MPIO_DEBUG(" Writing Pointdata... ")
784          // mark the active data arrays
785          bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
786          sprintf(header, "<PointData");
787          for (i_data =0 ;i_data<num_data;++i_data)
788          {
789            if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])
790            {
791              // if the rank == 0:   --> scalar data
792              // if the rank == 1:   --> vector data
793              // if the rank == 2:   --> tensor data
794    
795              switch(getDataPointRank(data_pp[i_data]))
796              {
797              case 0:
798                if (! set_scalar)
799                {
800                  sprintf(tmpBuf," Scalars=\"%s\"",names_p[i_data]);
801                  strcat(header,tmpBuf);
802                  set_scalar=TRUE;
803                }
804                break;
805              case 1:
806                if (! set_vector)
807                {
808                  sprintf(tmpBuf," Vectors=\"%s\"",names_p[i_data]);
809                  strcat(header,tmpBuf);
810                  set_vector=TRUE;
811                }
812                break;
813              case 2:
814                if (! set_tensor)
815                {
816                  sprintf(tmpBuf," Tensors=\"%s\"",names_p[i_data]);
817                  strcat(header,tmpBuf);
818                  set_tensor=TRUE;
819                }
820                break;
821              default:
822                sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);
823                Finley_setError(VALUE_ERROR,error_msg);
824                return;
825              }
826            }
827          }
828          strcat(header, ">\n");
829          MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);
830          MPI_Wait(&req,&status);
831        }
832      }
833    
834      // write actual data
835      if(write_pointdata)
836      {
837        for (i_data =0 ;i_data<num_data;++i_data)
838        {
839          if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])
840          {
841            numPointsPerSample = elements->ReferenceElement->numQuadNodes;
842            rank = getDataPointRank(data_pp[i_data]);
843            nComp = getDataPointSize(data_pp[i_data]);
844            nCompReqd=1;   // the number of components required by vtk
845            shape=0;
846            if (rank == 0)
847            {
848              nCompReqd = 1;
849            }
850            else if (rank == 1)
851            {
852              shape=getDataPointShape(data_pp[i_data], 0);
853              if  (shape>3)
854              {
855                Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");
856                return;
857              }
858              nCompReqd = 3;
859            }
860            else
861            {
862              shape=getDataPointShape(data_pp[i_data], 0);
863              if  (shape>3 || shape != getDataPointShape(data_pp[i_data], 1))
864              {
865                Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
866                return;
867              }
868              nCompReqd = 9;
869            }
870    
871            if( myRank == 0)
872            {
873              char header[250];
874              sprintf(header, "<DataArray Name=\"%s\" type=\"Float32\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);
875              MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);
876              MPI_Wait(&req,&status);
877            }
878            // write out the data
879            // if the number of required components is more than the number
880            // of actual components, pad with zeros
881    
882            char tmpbuf[14];
883            char* largebuf = MEMALLOC(nCompReqd*numLocalNodes*14 + numLocal*nCompReqd + nCompReqd + 14,char);
884            largebuf[0] = '\0';
885            bool_t do_write=TRUE;
886            size_t tsz = 0;
887    
888            for(k=0;k < nodeCache.size;k++)
889            {
890              i = nodeCache.values[k];
891    
892              if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
893              {
894                if (mesh_p->Nodes->toReduced[i]>=0)
895                {
896                  switch(getFunctionSpaceType(data_pp[i_data]))
897                  {
898                  case FINLEY_DEGREES_OF_FREEDOM:
899                    values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);
900                    break;
901                  case FINLEY_REDUCED_DEGREES_OF_FREEDOM:
902                    values = getSampleData(data_pp[i_data],mesh_p->Nodes->reducedDegreeOfFreedom[i]);
903                    break;
904                  case FINLEY_NODES:
905                    values = getSampleData(data_pp[i_data],i);
906                    break;
907                  }
908                  do_write=TRUE;
909                }
910                else
911                {
912                  do_write=FALSE;
913                }
914              }
915              else
916              {
917                do_write=TRUE;
918                switch(getFunctionSpaceType(data_pp[i_data]))
919                {
920                case FINLEY_DEGREES_OF_FREEDOM:
921                  values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);
922                  break;
923                case FINLEY_NODES:
924                  values = getSampleData(data_pp[i_data],i);
925                  break;
926                }
927              }
928              // write the data different ways for scalar, vector and tensor
929              if (do_write)
930              {
931                if (nCompReqd == 1)
932                {
933                  sprintf(tmpbuf," %e", values[0]);
934                  __len=strlen(tmpbuf);
935                  __STRCAT(largebuf,tmpbuf)
936                  tsz+=__len;
937                }
938                else if (nCompReqd == 3)
939                {
940                  for (m=0; m<shape; m++)
941                  {
942    
943                    sprintf(tmpbuf," %e",values[m]);
944                    __len=strlen(tmpbuf);
945                    __STRCAT(largebuf,tmpbuf)
946                    tsz+=__len;
947                  }
948                  for (m=0; m<nCompReqd-shape; m++)
949                  {
950                    __len=13;
951                    __STRCAT(largebuf,zero)
952                    tsz+=__len;
953                  }
954                }
955                else if (nCompReqd == 9)
956                {
957                  // tensor data, so have a 3x3 matrix to output as a row
958                  //  of 9 data points
959                  count = 0;
960                  for (m=0; m<shape; m++)
961                  {
962                    for (n=0; n<shape; n++)
963                    {
964                      sprintf(tmpbuf," %e", values[count]);
965                      __len=strlen(tmpbuf);
966                      __STRCAT(largebuf,tmpbuf)
967                      tsz+=__len;
968                      count++;
969                    }
970                    for (n=0; n<3-shape; n++)
971                    {
972                      __len=13;
973                      __STRCAT(largebuf,zero)
974                      tsz+=__len;
975                    }
976                  }
977                  for (m=0; m<3-shape; m++)
978                  {
979                    for (n=0; n<3; n++)
980                    {
981                      __len=13;
982                      __STRCAT(largebuf,zero)
983                      tsz+=__len;
984    
985                    }
986                  }
987                }
988                __len=1;
989                __STRCAT(largebuf,newline)
990                tsz+=1;
991              }
992    
993            }
994            // Write out local data
995            MPI_File_write_ordered(fh,largebuf,tsz,MPI_CHAR,&status);
996            MEMFREE(largebuf);
997            if( myRank == 0)
998            {
999              char *tag = "</DataArray>\n";
1000              MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);
1001              MPI_Wait(&req,&status);
1002            }
1003          }
1004        }
1005        // Finish off with closing tag
1006        if(myRank == 0)
1007        {
1008          char* tag =  "</PointData>\n";
1009          MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);
1010          MPI_Wait(&req,&status);
1011        }
1012        MPIO_DEBUG(" Done Writing Pointdata ")
1013      }
1014      // end write_pointdata
1015    
1016      // Write Cell data header Tags
1017      if(myRank == 0)
1018      {
1019        if( write_celldata)
1020        {
1021          char tmpBuf[80];
1022          char header[600];
1023          // mark the active data arrays
1024          bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
1025          sprintf(tmpBuf, "<CellData");
1026          strcat(header,tmpBuf);
1027          for (i_data =0 ;i_data<num_data;++i_data)
1028          {
1029            if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data])
1030            {
1031              // if the rank == 0:   --> scalar data
1032              // if the rank == 1:   --> vector data
1033              // if the rank == 2:   --> tensor data
1034    
1035              switch(getDataPointRank(data_pp[i_data]))
1036              {
1037              case 0:
1038                if (! set_scalar)
1039                {
1040                  sprintf(tmpBuf," Scalars=\"%s\"",names_p[i_data]);
1041                  strcat(header,tmpBuf);
1042                  set_scalar=TRUE;
1043                }
1044                break;
1045              case 1:
1046                if (! set_vector)
1047                {
1048                  sprintf(tmpBuf," Vectors=\"%s\"",names_p[i_data]);
1049                  set_vector=TRUE;
1050                }
1051                break;
1052              case 2:
1053                if (! set_tensor)
1054                {
1055                  sprintf(tmpBuf," Tensors=\"%s\"",names_p[i_data]);
1056                  set_tensor=TRUE;
1057                }
1058                break;
1059              default:
1060                sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);
1061                Finley_setError(VALUE_ERROR,error_msg);
1062                return;
1063              }
1064            }
1065          }
1066          strcat(header, ">\n");
1067          MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);
1068          MPI_Wait(&req,&status);
1069        }
1070      }
1071    
1072      // write actual data (collective)
1073      if(write_celldata)
1074      {
1075        for (i_data =0 ;i_data<num_data;++i_data)
1076        {
1077          if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data])
1078          {
1079            numPointsPerSample = elements->ReferenceElement->numQuadNodes;
1080            rank = getDataPointRank(data_pp[i_data]);
1081            nComp = getDataPointSize(data_pp[i_data]);
1082            nCompReqd=1;   // the number of components required by vtk
1083            shape=0;
1084            if (rank == 0)
1085            {
1086              nCompReqd = 1;
1087            }
1088            else if (rank == 1)
1089            {
1090              shape=getDataPointShape(data_pp[i_data], 0);
1091              if  (shape>3)
1092              {
1093                Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");
1094                return;
1095              }
1096              nCompReqd = 3;
1097            }
1098            else
1099            {
1100              shape=getDataPointShape(data_pp[i_data], 0);
1101              if  (shape>3 || shape != getDataPointShape(data_pp[i_data], 1))
1102              {
1103                Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
1104                return;
1105              }
1106              nCompReqd = 9;
1107            }
1108    
1109            if( myRank == 0)
1110            {
1111              char header[250];
1112              sprintf(header,"<DataArray Name=\"%s\" type=\"Float32\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);
1113              MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);
1114              MPI_Wait(&req,&status);
1115            }
1116    
1117            // Write the actual data */
1118            char tmpbuf[14];
1119            char* largebuf = MEMALLOC(nCompReqd*numLocalCells*14 + numLocalCells*nCompReqd + nCompReqd + 14,char);
1120            largebuf[0] = '\0';
1121            size_t tsz = 0;
1122    
1123            double sampleAvg[nComp];
1124    
1125            for (k=0; i<elementCache.size; k++)
1126            {
1127              i = elementCache.values[k];
1128    
1129              values = getSampleData(data_pp[i_data], i);
1130              // averaging over the number of points in the sample
1131              for (k=0; k<nComp; k++)
1132              {
1133                rtmp = 0.;
1134                for (j=0; j<numPointsPerSample; j++) rtmp += values[INDEX2(k,j,nComp)];
1135                sampleAvg[k] = rtmp/numPointsPerSample;
1136              }
1137              // if the number of required components is more than the number
1138              // of actual components, pad with zeros
1139    
1140              // probably only need to get shape of first element
1141              // write the data different ways for scalar, vector and tensor
1142              if (nCompReqd == 1)
1143              {
1144                sprintf(tmpbuf, " %e", sampleAvg[0]);
1145                __len=strlen(tmpbuf);
1146                __STRCAT(largebuf,tmpbuf)
1147                tsz+=__len;
1148              }
1149              else if (nCompReqd == 3)
1150              {
1151                // write out the data
1152                for (m=0; m<shape; m++)
1153                {
1154                  sprintf(tmpbuf, " %e", sampleAvg[m]);
1155                  __len=strlen(tmpbuf);
1156                  __STRCAT(largebuf,tmpbuf)
1157                  tsz+=__len;
1158                }
1159                for (m=0; m<nCompReqd-shape; m++)
1160                {
1161                  tsz+=13;
1162                  strcat(largebuf," 0.000000e+00");
1163    
1164                }
1165              }
1166              else if (nCompReqd == 9)
1167              {
1168                // tensor data, so have a 3x3 matrix to output as a row
1169                // of 9 data points
1170                count = 0;
1171                for (m=0; m<shape; m++)
1172                {
1173                  for (n=0; n<shape; n++)
1174                  {
1175                    sprintf(tmpbuf, " %e", sampleAvg[count]);
1176                    __len=strlen(tmpbuf);
1177                    __STRCAT(largebuf,tmpbuf)
1178                    tsz+=__len;;
1179    
1180    
1181                    count++;
1182                  }
1183                  for (n=0; n<3-shape; n++)
1184                  {
1185                    __len=13;
1186                    __STRCAT(largebuf,zero)
1187                    tsz+=13;
1188                  }
1189                }
1190                for (m=0; m<3-shape; m++)
1191                  for (n=0; n<3; n++)
1192                  {
1193                    __len=13;
1194                    __STRCAT(largebuf,zero)
1195                    tsz+=13;
1196                  }
1197              }
1198              __len=1;
1199              __STRCAT(largebuf,newline)
1200              tsz+=1;
1201            }
1202            MPI_File_write_ordered(fh,largebuf,tsz,MPI_CHAR,&status);
1203            MEMFREE(largebuf);
1204            if( myRank == 0)
1205            {
1206              char *tag = "</DataArray>\n";
1207              MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);
1208              MPI_Wait(&req,&status);
1209            }
1210    
1211          }
1212        }
1213        // closing celldata tag
1214        if(myRank == 0)
1215        {
1216          char* tag =  "</CellData>\n";
1217          MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);
1218          MPI_Wait(&req,&status);
1219        }
1220      }
1221    
1222      /* tag and bag... */
1223      if (myRank == 0)
1224      {
1225        char *footer = "</Piece>\n</UnstructuredGrid>\n</VTKFile>";
1226        MPI_File_iwrite_shared(fh,footer,strlen(footer),MPI_CHAR,&req);
1227        MPI_Wait(&req,&status);
1228      }
1229    
1230      MEMFREE(nodesGlobal);
1231      MEMFREE(nodeCache.values);
1232      MEMFREE(elementCache.values);
1233    #ifdef MPIO_HINTS
1234      MPI_Info_free(&infoHints);
1235    #undef MPIO_HINTS
1236    #endif
1237      printf("\nTime: %f \n",  Paso_timer() - time0);
1238      MPI_File_close(&fh);
1239      MPIO_DEBUG(" ***** Exit saveVTK ***** ")
1240    #undef __STRCAT
1241    }
1242    
1243    #undef MPIO_DEBUG
1244    #else
1245    
1246    
1247    
1248    
1249    void Finley_Mesh_saveVTK(const char * filename_p, Finley_Mesh *mesh_p, const dim_t num_data,char* *names_p, escriptDataC* *data_pp)
1250    {
1251      char error_msg[LenErrorMsg_MAX];
1252      /* if there is no mesh we just return */
1253      if (mesh_p==NULL) return;
1254    
1255      int i, j, k, numVTKNodesPerElement,i_data,m, count, n, rank,shape, numPoints, cellType, numCells,
1256      nDim, numPointsPerSample, nComp, nCompReqd;
1257    
1258      index_t j2;
1259      double* values, rtmp;
1260      char elemTypeStr[32];
1261    
1262    /* open the file and check handle */    /* open the file and check handle */
1263    
1264    FILE * fileHandle_p = fopen(filename_p, "w");    FILE * fileHandle_p = fopen(filename_p, "w");
1265    if (fileHandle_p==NULL) {    if (fileHandle_p==NULL)
1266      Finley_ErrorCode=IO_ERROR;    {
1267      sprintf(Finley_ErrorMsg,      sprintf(error_msg, "saveVTK: File %s could not be opened for writing.", filename_p);
1268          "File %s could not be opened for writing.", filename_p);      Finley_setError(IO_ERROR,error_msg);
1269        return;
1270      }
1271      /* find the mesh type to be written */
1272      int nodetype=FINLEY_DEGREES_OF_FREEDOM;
1273      int elementtype=FINLEY_UNKNOWN;
1274      bool_t isCellCentered[num_data],write_celldata=FALSE,write_pointdata=FALSE;
1275      for (i_data=0;i_data<num_data;++i_data)
1276      {
1277        if (! isEmpty(data_pp[i_data]))
1278        {
1279          switch(getFunctionSpaceType(data_pp[i_data]))
1280          {
1281          case FINLEY_DEGREES_OF_FREEDOM:
1282            nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1283            if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)
1284            {
1285              elementtype=FINLEY_ELEMENTS;
1286            }
1287            else
1288            {
1289              Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1290              fclose(fileHandle_p);
1291              return;
1292            }
1293            isCellCentered[i_data]=FALSE;
1294            break;
1295          case FINLEY_REDUCED_DEGREES_OF_FREEDOM:
1296            nodetype = FINLEY_REDUCED_DEGREES_OF_FREEDOM;
1297            if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)
1298            {
1299              elementtype=FINLEY_ELEMENTS;
1300            }
1301            else
1302            {
1303              Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1304              fclose(fileHandle_p);
1305              return;
1306            }
1307            isCellCentered[i_data]=FALSE;
1308            break;
1309          case FINLEY_NODES:
1310            nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1311            if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)
1312            {
1313              elementtype=FINLEY_ELEMENTS;
1314            }
1315            else
1316            {
1317              Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1318              fclose(fileHandle_p);
1319              return;
1320            }
1321            isCellCentered[i_data]=FALSE;
1322            break;
1323          case FINLEY_ELEMENTS:
1324            nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1325            if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)
1326            {
1327              elementtype=FINLEY_ELEMENTS;
1328            }
1329            else
1330            {
1331              Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1332              fclose(fileHandle_p);
1333              return;
1334            }
1335            isCellCentered[i_data]=TRUE;
1336            break;
1337          case FINLEY_FACE_ELEMENTS:
1338            nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1339            if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_FACE_ELEMENTS)
1340            {
1341              elementtype=FINLEY_FACE_ELEMENTS;
1342            }
1343            else
1344            {
1345              Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1346              fclose(fileHandle_p);
1347              return;
1348            }
1349            isCellCentered[i_data]=TRUE;
1350            break;
1351          case FINLEY_POINTS:
1352            nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1353            if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_POINTS)
1354            {
1355              elementtype=FINLEY_POINTS;
1356            }
1357            else
1358            {
1359              Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1360              fclose(fileHandle_p);
1361              return;
1362            }
1363            isCellCentered[i_data]=TRUE;
1364            break;
1365          case FINLEY_CONTACT_ELEMENTS_1:
1366            nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1367            if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1)
1368            {
1369              elementtype=FINLEY_CONTACT_ELEMENTS_1;
1370            }
1371            else
1372            {
1373              Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1374              fclose(fileHandle_p);
1375              return;
1376            }
1377            isCellCentered[i_data]=TRUE;
1378            break;
1379          case FINLEY_CONTACT_ELEMENTS_2:
1380            nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1381            if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1)
1382            {
1383              elementtype=FINLEY_CONTACT_ELEMENTS_1;
1384            }
1385            else
1386            {
1387              Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1388              fclose(fileHandle_p);
1389              return;
1390            }
1391            isCellCentered[i_data]=TRUE;
1392            break;
1393          default:
1394            sprintf(error_msg,"saveVTK: Finley does not know anything about function space type %d",getFunctionSpaceType(data_pp[i_data]));
1395            Finley_setError(TYPE_ERROR,error_msg);
1396            fclose(fileHandle_p);
1397            return;
1398          }
1399          if (isCellCentered[i_data])
1400          {
1401            write_celldata=TRUE;
1402          }
1403          else
1404          {
1405            write_pointdata=TRUE;
1406          }
1407        }
1408      }
1409      /* select nomber of points and the mesh component */
1410      numPoints = mesh_p->Nodes->numNodes;
1411      if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1412      {
1413        numPoints = mesh_p->Nodes->reducedNumNodes;
1414      }
1415      else
1416      {
1417        numPoints = mesh_p->Nodes->numNodes;
1418      }
1419      if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;
1420      Finley_ElementFile* elements=NULL;
1421      switch(elementtype)
1422      {
1423      case FINLEY_ELEMENTS:
1424        elements=mesh_p->Elements;
1425        break;
1426      case FINLEY_FACE_ELEMENTS:
1427        elements=mesh_p->FaceElements;
1428        break;
1429      case FINLEY_POINTS:
1430        elements=mesh_p->Points;
1431        break;
1432      case FINLEY_CONTACT_ELEMENTS_1:
1433        elements=mesh_p->ContactElements;
1434        break;
1435      }
1436      if (elements==NULL)
1437      {
1438        Finley_setError(SYSTEM_ERROR,"saveVTK: undefined element file");
1439        fclose(fileHandle_p);
1440        return;
1441      }
1442      /* map finley element type to VTK element type */
1443      numCells = elements->numElements;
1444      ElementTypeId TypeId;
1445      if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1446      {
1447        TypeId = elements->LinearReferenceElement->Type->TypeId;
1448      }
1449      else
1450      {
1451        TypeId = elements->ReferenceElement->Type->TypeId;
1452      }
1453    
1454      switch(TypeId)
1455      {
1456      case Point1:
1457      case Line2Face:
1458      case Line3Face:
1459      case Point1_Contact:
1460      case Line2Face_Contact:
1461      case Line3Face_Contact:
1462        cellType = VTK_VERTEX;
1463        numVTKNodesPerElement = 1;
1464        strcpy(elemTypeStr, "VTK_VERTEX");
1465        break;
1466    
1467      case Line2:
1468      case Tri3Face:
1469      case Rec4Face:
1470      case Line2_Contact:
1471      case Tri3_Contact:
1472      case Tri3Face_Contact:
1473      case Rec4Face_Contact:
1474        cellType = VTK_LINE;
1475        numVTKNodesPerElement = 2;
1476        strcpy(elemTypeStr, "VTK_LINE");
1477        break;
1478    
1479      case Tri3:
1480      case Tet4Face:
1481      case Tet4Face_Contact:
1482        cellType = VTK_TRIANGLE;
1483        numVTKNodesPerElement = 3;
1484        strcpy(elemTypeStr, "VTK_TRIANGLE");
1485        break;
1486    
1487      case Rec4:
1488      case Hex8Face:
1489      case Rec4_Contact:
1490      case Hex8Face_Contact:
1491        cellType = VTK_QUAD;
1492        numVTKNodesPerElement = 4;
1493        strcpy(elemTypeStr, "VTK_QUAD");
1494        break;
1495    
1496      case Tet4:
1497        cellType = VTK_TETRA;
1498        numVTKNodesPerElement = 4;
1499        strcpy(elemTypeStr, "VTK_TETRA");
1500        break;
1501    
1502      case Hex8:
1503        cellType = VTK_HEXAHEDRON;
1504        numVTKNodesPerElement = 8;
1505        strcpy(elemTypeStr, "VTK_HEXAHEDRON");
1506        break;
1507    
1508      case Line3:
1509      case Tri6Face:
1510      case Rec8Face:
1511      case Line3_Contact:
1512      case Tri6Face_Contact:
1513      case Rec8Face_Contact:
1514        cellType = VTK_QUADRATIC_EDGE;
1515        numVTKNodesPerElement = 3;
1516        strcpy(elemTypeStr, "VTK_QUADRATIC_EDGE");
1517        break;
1518    
1519      case Tri6:
1520      case Tet10Face:
1521      case Tri6_Contact:
1522      case Tet10Face_Contact:
1523        cellType = VTK_QUADRATIC_TRIANGLE;
1524        numVTKNodesPerElement = 6;
1525        strcpy(elemTypeStr, "VTK_QUADRATIC_TRIANGLE");
1526        break;
1527    
1528      case Rec8:
1529      case Hex20Face:
1530      case Rec8_Contact:
1531      case Hex20Face_Contact:
1532        cellType = VTK_QUADRATIC_QUAD;
1533        numVTKNodesPerElement = 8;
1534        strcpy(elemTypeStr, "VTK_QUADRATIC_QUAD");
1535        break;
1536    
1537      case Tet10:
1538        cellType = VTK_QUADRATIC_TETRA;
1539        numVTKNodesPerElement = 10;
1540        strcpy(elemTypeStr, "VTK_QUADRATIC_TETRA");
1541        break;
1542    
1543      case Hex20:
1544        cellType = VTK_QUADRATIC_HEXAHEDRON;
1545        numVTKNodesPerElement = 20;
1546        strcpy(elemTypeStr, "VTK_QUADRATIC_HEXAHEDRON");
1547        break;
1548    
1549      default:
1550        sprintf(error_msg, "saveVTK: Element type %s is not supported by VTK",elements->ReferenceElement->Type->Name);
1551        Finley_setError(VALUE_ERROR,error_msg);
1552        fclose(fileHandle_p);
1553      return;      return;
1554    }    }
1555    /* xml header */    /* xml header */
1556    fprintf(fileHandle_p, "<?xml version=\"1.0\"?>\n");    fprintf(fileHandle_p, "<?xml version=\"1.0\"?>\n");
1557    fprintf(fileHandle_p,    fprintf(fileHandle_p, "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\">\n");
       "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\">\n");  
1558    
1559    /* finley uses an unstructured mesh, so UnstructuredGrid *should* work */    /* finley uses an unstructured mesh, so UnstructuredGrid *should* work */
1560    fprintf(fileHandle_p, "<UnstructuredGrid>\n");    fprintf(fileHandle_p, "<UnstructuredGrid>\n");
1561    
1562    /* is there only one "piece" to the data?? */    /* is there only one "piece" to the data?? */
1563    fprintf(fileHandle_p, "<Piece "    fprintf(fileHandle_p, "<Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n",numPoints, numCells);
       "NumberOfPoints=\"%d\" "  
       "NumberOfCells=\"%d\">\n",  
       numPoints, numCells);  
   
1564    /* now for the points; equivalent to positions section in saveDX() */    /* now for the points; equivalent to positions section in saveDX() */
1565    /* "The points element explicitly defines coordinates for each point    /* "The points element explicitly defines coordinates for each point
1566     * individually.  It contains one DataArray element describing an array    * individually.  It contains one DataArray element describing an array
1567     * with three components per value, each specifying the coordinates of one    * with three components per value, each specifying the coordinates of one
1568     * point" - from Vtk User's Guide    * point" - from Vtk User's Guide
1569     */    */
1570    fprintf(fileHandle_p, "<Points>\n");    fprintf(fileHandle_p, "<Points>\n");
1571    /*    /*
1572     * the reason for this if statement is explained in the long comment below    * the reason for this if statement is explained in the long comment below
1573     */    */
1574    if (nDim < 3) {    nDim = mesh_p->Nodes->numDim;
1575      fprintf(fileHandle_p, "<DataArray "    fprintf(fileHandle_p, "<DataArray NumberOfComponents=\"%d\" type=\"Float32\" format=\"ascii\">\n",MAX(3,nDim));
1576          "NumberOfComponents=\"3\" "    /* vtk/mayavi doesn't like 2D data, it likes 3D data with a degenerate
1577          "type=\"Float32\" "    * third dimension to handle 2D data (like a sheet of paper).  So, if
1578          "format=\"ascii\">\n");    * nDim is 2, we have to append zeros to the array to get this third
1579    } else {    * dimension, and keep the visualisers happy.
1580      fprintf(fileHandle_p, "<DataArray "    * Indeed, if nDim is less than 3, must pad all empty dimensions, so
1581          "NumberOfComponents=\"%d\" "    * that the total number of dims is 3.
1582          "type=\"Float32\" "    */
1583          "format=\"ascii\">\n",    if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1584          nDim);    {
1585    }      for (i = 0; i < mesh_p->Nodes->numNodes; i++)
1586    for (i = 0; i < mesh_p->Nodes->numNodes; i++) {      {
1587      fprintf(fileHandle_p,        if (mesh_p->Nodes->toReduced[i]>=0)
1588          "%e", mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)]);        {
1589      for (j = 1; j < nDim; j++) {          for (j = 0; j < nDim; j++) fprintf(fileHandle_p, " %e",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);
1590        fprintf(fileHandle_p,          for (k=0; k<3-nDim; k++) fprintf(fileHandle_p, " %e",0.);
1591            " %f",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);          fprintf(fileHandle_p, "\n");
       /* 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 (nDim < 3) {  
     for (k=0; k<3-nDim; k++) {  
       fprintf(fileHandle_p, " 0");  
     }  
1592        }        }
1593      }      }
1594      fprintf(fileHandle_p, "\n");    }
1595    }    else
1596      {
1597        for (i = 0; i < mesh_p->Nodes->numNodes; i++)
1598        {
1599    
1600          for (j = 0; j < nDim; j++) fprintf(fileHandle_p, " %e",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);
1601          for (k=0; k<3-nDim; k++) fprintf(fileHandle_p, " %e",0.);
1602          fprintf(fileHandle_p, "\n");
1603        }
1604    
1605      }
1606    fprintf(fileHandle_p, "</DataArray>\n");    fprintf(fileHandle_p, "</DataArray>\n");
1607    fprintf(fileHandle_p, "</Points>\n");    fprintf(fileHandle_p, "</Points>\n");
1608    
1609    /* connections */    /* write out the DataArray element for the connectivity */
1610    /* now for the cells */  
1611    /* "The Cells element defines cells explicitly by specifying point    int NN = elements->ReferenceElement->Type->numNodes;
1612     * connectivity and cell types.  It contains three DataArray elements.  The    fprintf(fileHandle_p, "<Cells>\n");
1613     * first array specifies the point connectivity.  All cells' point lists    fprintf(fileHandle_p, "<DataArray Name=\"connectivity\" type=\"Int32\" format=\"ascii\">\n");
1614     * are concatenated together.  The second array specifies th eoffset into  
1615     * the connectivity array for the end of each cell.  The third array    if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1616     * specifies the type of each cell.    {
1617     */      for (i = 0; i < numCells; i++)
1618    /* if no element table is present jump over the connection table */      {
1619    int cellType;        for (j = 0; j < numVTKNodesPerElement; j++)
1620    if (elements!=NULL) {          fprintf(fileHandle_p,"%d ",mesh_p->Nodes->toReduced[elements->Nodes[INDEX2(elements->ReferenceElement->Type->linearNodes[j], i, NN)]]);
1621      fprintf(fileHandle_p, "<Cells>\n");        fprintf(fileHandle_p, "\n");
     ElementTypeId TypeId = elements->ReferenceElement->Type->TypeId;  
     switch(TypeId) {  
     case Point1:  
       cellType = VTK_VERTEX;  
       break;  
     case Line2:  
       cellType = VTK_LINE;  
       break;  
     case Line3:  
       cellType = VTK_QUADRATIC_EDGE;  
       break;  
     case Tri3:  
       cellType = VTK_TRIANGLE;  
       break;  
     case Tri6:  
       cellType = VTK_QUADRATIC_TRIANGLE;  
       break;  
     case Rec4:  
       cellType = VTK_QUAD;  
       break;  
     case Rec8:  
       cellType = VTK_QUADRATIC_QUAD;  
       break;  
     case Tet4:  
       cellType = VTK_TETRA;  
       break;  
     case Tet10:  
       cellType = VTK_QUADRATIC_TETRA;  
       break;  
     case Hex8:  
       cellType = VTK_HEXAHEDRON;  
       break;  
     case Hex20:  
       cellType = VTK_QUADRATIC_HEXAHEDRON;  
       break;  
     case Line2Face:  
       cellType = VTK_VERTEX;  
       break;  
     case Line3Face:  
       cellType = VTK_VERTEX;  
       break;  
     case Tri3Face:  
       cellType = VTK_LINE;  
       break;  
     case Tri6Face:  
       cellType = VTK_QUADRATIC_EDGE;  
       break;  
     case Rec4Face:  
       cellType = VTK_LINE;  
       break;  
     case Rec8Face:  
       cellType = VTK_QUADRATIC_EDGE;  
       break;  
     case Tet4Face:  
       cellType = VTK_TRIANGLE;  
       break;  
     case Tet10Face:  
       cellType = VTK_QUADRATIC_TRIANGLE;  
       break;  
     case Hex8Face:  
       cellType = VTK_QUADRATIC_QUAD;  
       break;  
     case Hex20Face:  
       cellType = VTK_QUADRATIC_QUAD;  
       break;  
     case Point1_Contact:  
       cellType = VTK_VERTEX;  
       break;  
     case Line2_Contact:  
       cellType = VTK_LINE;  
       break;  
     case Line3_Contact:  
       cellType = VTK_QUADRATIC_EDGE;  
       break;  
     case Tri3_Contact:  
       cellType = VTK_TRIANGLE;  
       break;  
     case Tri6_Contact:  
       cellType = VTK_QUADRATIC_TRIANGLE;  
       break;  
     case Rec4_Contact:  
       cellType = VTK_QUAD;  
       break;  
     case Rec8_Contact:  
       cellType = VTK_QUADRATIC_QUAD;  
       break;  
     case Line2Face_Contact:  
       cellType = VTK_VERTEX;  
       break;  
     case Line3Face_Contact:  
       cellType = VTK_VERTEX;  
       break;  
     case Tri3Face_Contact:  
       cellType = VTK_LINE;  
       break;  
     case Tri6Face_Contact:  
       cellType = VTK_QUADRATIC_EDGE;  
       break;  
     case Rec4Face_Contact:  
       cellType = VTK_LINE;  
       break;  
     case Rec8Face_Contact:  
       cellType = VTK_QUADRATIC_EDGE;  
       break;  
     case Tet4Face_Contact:  
       cellType = VTK_TRIANGLE;  
       break;  
     case Tet10Face_Contact:  
       cellType = VTK_QUADRATIC_TRIANGLE;  
       break;  
     case Hex8Face_Contact:  
       cellType = VTK_QUAD;  
       break;  
     case Hex20Face_Contact:  
       cellType = VTK_QUADRATIC_QUAD;  
       break;  
     default:  
       Finley_ErrorCode=VALUE_ERROR;  
       sprintf(Finley_ErrorMsg,  
           "Element type %s is not supported by VTK",  
           elements->ReferenceElement->Type->Name);  
       return;  
     }  
   
     switch(cellType) {  
     case VTK_VERTEX:  
       numVTKNodesPerElement = 1;  
       strcpy(elemTypeStr, "VTK_VERTEX");  
       break;  
     case VTK_LINE:  
       numVTKNodesPerElement = 2;  
       strcpy(elemTypeStr, "VTK_LINE");  
       break;  
     case VTK_TRIANGLE:  
       numVTKNodesPerElement = 3;  
       strcpy(elemTypeStr, "VTK_TRIANGLE");  
       break;  
     case VTK_QUAD:  
       numVTKNodesPerElement = 4;  
       strcpy(elemTypeStr, "VTK_QUAD");  
       break;  
     case VTK_TETRA:  
       numVTKNodesPerElement = 4;  
       strcpy(elemTypeStr, "VTK_TETRA");  
       break;  
     case VTK_HEXAHEDRON:  
       numVTKNodesPerElement = 8;  
       strcpy(elemTypeStr, "VTK_HEXAHEDRON");  
       break;  
     case VTK_QUADRATIC_EDGE:  
       numVTKNodesPerElement = 3;  
       strcpy(elemTypeStr, "VTK_QUADRATIC_EDGE");  
       break;  
     case VTK_QUADRATIC_TRIANGLE:  
       numVTKNodesPerElement = 6;  
       strcpy(elemTypeStr, "VTK_QUADRATIC_TRIANGLE");  
       break;  
     case VTK_QUADRATIC_QUAD:  
       numVTKNodesPerElement = 8;  
       strcpy(elemTypeStr, "VTK_QUADRATIC_QUAD");  
       break;  
     case VTK_QUADRATIC_TETRA:  
       numVTKNodesPerElement = 10;  
       strcpy(elemTypeStr, "VTK_QUADRATIC_TETRA");  
       break;  
     case VTK_QUADRATIC_HEXAHEDRON:  
       numVTKNodesPerElement = 20;  
       strcpy(elemTypeStr, "VTK_QUADRATIC_HEXAHEDRON");  
       break;  
     default:  
       Finley_ErrorCode = VALUE_ERROR;  
       sprintf(Finley_ErrorMsg,  
           "Cell type %d is not supported by VTK", cellType);  
       return;  
     }  
   
     /* write out the DataArray element for the connectivity */  
     fprintf(fileHandle_p, "<DataArray "  
         "Name=\"connectivity\" "  
         "type=\"Int32\" "  
         "format=\"ascii\">\n");  
     int NN = elements->ReferenceElement->Type->numNodes;  
     int counter = 0;  
     for (i = 0; i < numCells; i++) {  
       fprintf(fileHandle_p, "%d ",  
           mesh_p->Elements->Nodes[INDEX2(0, i, NN)]);  
       counter++; /* counter for the number of connectivity points written */  
       /* if the counter gets too big (i.e. the line gets too long),  
        * then add a newline and reset */  
       if (counter > 19) {  
       fprintf(fileHandle_p, "\n");  
       counter = 0;  
       }  
       for (j = 1; j < numVTKNodesPerElement; j++) {  
     fprintf(fileHandle_p,"%d ", mesh_p->Elements->Nodes[INDEX2(j, i, NN)]);  
     counter++;  
     /* if the counter gets too big (i.e. the line gets too long),  
      * then add a newline and reset */  
     if (counter > 19) {  
         fprintf(fileHandle_p, "\n");  
         counter = 0;  
     }  
       }  
     }  
     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");  
     counter = 0;  /* counter for the number of offsets written to file */  
     for (i=numVTKNodesPerElement; i<=numCells*numVTKNodesPerElement; i+=numVTKNodesPerElement) {  
       fprintf(fileHandle_p, "%d ", i);  
       counter++;  
       /* if the counter gets too big (i.e. the line gets too long),  
        * then add a newline and reset */  
       if (counter > 19) {  
       counter = 0;  
       fprintf(fileHandle_p, "\n");  
       }  
     }  
     fprintf(fileHandle_p, "\n");  
     fprintf(fileHandle_p, "</DataArray>\n");  
   
     /* write out the DataArray element for the types */  
     counter = 0; /* counter for the number of types written to file */  
     fprintf(fileHandle_p, "<DataArray "  
         "Name=\"types\" "  
         "type=\"UInt8\" "  
         "format=\"ascii\">\n");  
     for (i=0; i<numCells; i++) {  
       fprintf(fileHandle_p, "%d ", cellType);  
       counter++;  
       /* if the counter gets too big (i.e. the line gets too long),  
        * then add a newline and reset */  
       if (counter > 30) {  
       counter = 0;  
       fprintf(fileHandle_p, "\n");  
       }  
     }  
     fprintf(fileHandle_p, "\n");  
     fprintf(fileHandle_p, "</DataArray>\n");  
   
     /* finish off the <Cells> element */  
     fprintf(fileHandle_p, "</Cells>\n");  
   }  
   
   /* data */  
   if (!isEmpty(data_p)) {  
     int rank = getDataPointRank(data_p);  
     int nComp = getDataPointSize(data_p);  
     /* barf if rank is greater than two */  
     if (rank > 2) {  
       Finley_ErrorCode = VALUE_ERROR;  
       sprintf(Finley_ErrorMsg,  
           "Vtk can't handle objects with rank greater than 2\n"  
           "object rank = %d\n", rank);  
       return;  
     }  
     /* if the rank == 0:   --> scalar data  
      * if the rank == 1:   --> vector data  
      * if the rank == 2:   --> tensor data  
      */  
     char dataNameStr[31], dataTypeStr[63];  
     int nCompReqd;   /* the number of components required by vtk */  
     if (rank == 0) {  
       strcpy(dataNameStr, "escript_scalar_data");  
       sprintf(dataTypeStr, "Scalars=\"%s\"", dataNameStr);  
       nCompReqd = 1;  
     }  
     else if (rank == 1) {  
       strcpy(dataNameStr, "escript_vector_data");  
       sprintf(dataTypeStr, "Vectors=\"%s\"", dataNameStr);  
       nCompReqd = 3;  
     }  
     else if (rank == 2) {  
       strcpy(dataNameStr, "escript_tensor_data");  
       sprintf(dataTypeStr, "Tensors=\"%s\"", dataNameStr);  
       nCompReqd = 9;  
     }  
     /* if have cell centred data then use CellData element,  
      * if have node centred data, then use PointData element  
      */  
     if (isCellCentered) {  
       /* now for the cell data */  
       fprintf(fileHandle_p, "<CellData %s>\n", dataTypeStr);  
       fprintf(fileHandle_p,  
           "<DataArray "  
           "Name=\"%s\" "  
           "type=\"Float32\" "  
           "NumberOfComponents=\"%d\" "  
           "format=\"ascii\">\n",  
           dataNameStr, nCompReqd);  
       int numPointsPerSample = elements->ReferenceElement->numQuadNodes;  
       if (numPointsPerSample) {  
     for (i=0; i<numCells; i++) {  
       values = getSampleData(data_p, i);  
       double sampleAvg[nComp];  
       for (k=0; k<nComp; k++) {  
         /* averaging over the number of points in the sample */  
         rtmp = 0.;  
         for (j=0; j<numPointsPerSample; j++) {  
           rtmp += values[INDEX2(k,j,nComp)];  
         }  
         sampleAvg[k] = rtmp/numPointsPerSample;  
       }  
       /* if the number of required components is more than the number  
        * of actual components, pad with zeros  
        */  
       /* probably only need to get shape of first element */  
       int shape = getDataPointShape(data_p, 0);  
       if (shape > 3) {  
         Finley_ErrorCode = VALUE_ERROR;  
         sprintf(Finley_ErrorMsg,  
             "shape should be 1, 2, or 3; I got %d\n", shape);  
         return;  
       }  
       /* write the data different ways for scalar, vector and tensor */  
       if (nCompReqd == 1) {  
         fprintf(fileHandle_p, " %f", sampleAvg[0]);  
       }  
       else if (nCompReqd == 3) {  
         int shape = getDataPointShape(data_p, 0);  
         /* write out the data */  
         for (int m=0; m<shape; m++) {  
           fprintf(fileHandle_p, " %f", sampleAvg[m]);  
         }  
         /* pad with zeros */  
         for (int m=0; m<nCompReqd-shape; m++) {  
           fprintf(fileHandle_p, " 0");  
         }  
       }  
       else if (nCompReqd == 9) {  
         /* tensor data, so have a 3x3 matrix to output as a row  
          * of 9 data points */  
         int count = 0;  
         int elems = 0;  
         for (int m=0; m<shape; m++) {  
           for (int n=0; n<shape; n++) {  
         fprintf(fileHandle_p, " %f", sampleAvg[count]);  
         count++;  
         elems++;  
           }  
           for (int n=0; n<3-shape; n++) {  
         fprintf(fileHandle_p, " 0");  
         elems++;  
           }  
         }  
         for (int m=0; m<nCompReqd-elems; m++) {  
           fprintf(fileHandle_p, " 0");  
         }  
       }  
       fprintf(fileHandle_p, "\n");  
     }  
       }  
       fprintf(fileHandle_p, "</DataArray>\n");  
       fprintf(fileHandle_p, "</CellData>\n");  
     } else {  
       /* now for the point data */  
       fprintf(fileHandle_p, "<PointData %s>\n", dataTypeStr);  
       fprintf(fileHandle_p, "<DataArray "  
           "Name=\"%s\" "  
           "type=\"Float32\" "  
           "NumberOfComponents=\"%d\" "  
           "format=\"ascii\">\n",  
           dataNameStr, nCompReqd);  
       for (i=0; i<mesh_p->Nodes->numNodes; i++) {  
     switch (nodetype) {  
     case(FINLEY_DEGREES_OF_FREEDOM):  
       values = getSampleData(data_p,  
                  mesh_p->Nodes->degreeOfFreedom[i]);  
       break;  
     case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):  
       if (mesh_p->Nodes->toReduced[i]>=0) {  
         values = getSampleData(data_p,  
                    mesh_p->Nodes->reducedDegreeOfFreedom[i]);  
       }  
       break;  
     case(FINLEY_NODES):  
       values = getSampleData(data_p,i);  
       break;  
     }  
     /* write out the data */  
     /* if the number of required components is more than the number  
      * of actual components, pad with zeros  
      */  
     /* probably only need to get shape of first element */  
     int shape = getDataPointShape(data_p, 0);  
     if (shape > 3) {  
       Finley_ErrorCode = VALUE_ERROR;  
       sprintf(Finley_ErrorMsg,  
           "shape should be 1, 2, or 3; I got %d\n", shape);  
       return;  
     }  
     /* write the data different ways for scalar, vector and tensor */  
     if (nCompReqd == 1) {  
       fprintf(fileHandle_p, " %f", values[0]);  
     }  
     else if (nCompReqd == 3) {  
       int shape = getDataPointShape(data_p, 0);  
       /* write out the data */  
       for (int m=0; m<shape; m++) {  
         fprintf(fileHandle_p, " %f", values[m]);  
       }  
       /* pad with zeros */  
       for (int m=0; m<nCompReqd-shape; m++) {  
         fprintf(fileHandle_p, " 0");  
       }  
     }  
     else if (nCompReqd == 9) {  
       /* tensor data, so have a 3x3 matrix to output as a row  
        * of 9 data points */  
       int count = 0;  
       int elems = 0;  
       for (int m=0; m<shape; m++) {  
         for (int n=0; n<shape; n++) {  
           fprintf(fileHandle_p, " %f", values[count]);  
           count++;  
           elems++;  
         }  
         for (int n=0; n<3-shape; n++) {  
           fprintf(fileHandle_p, " 0");  
           elems++;  
         }  
       }  
       for (int m=0; m<nCompReqd-elems; m++) {  
         fprintf(fileHandle_p, " 0");  
       }  
     }  
     fprintf(fileHandle_p, "\n");  
       }  
       fprintf(fileHandle_p, "</DataArray>\n");  
       fprintf(fileHandle_p, "</PointData>\n");  
1622      }      }
1623    }    }
1624      else if (VTK_QUADRATIC_HEXAHEDRON==cellType)
1625      {
1626        for (i = 0; i < numCells; i++)
1627        {
1628          fprintf(fileHandle_p,"%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d\n",
1629                  elements->Nodes[INDEX2(0, i, NN)],
1630                  elements->Nodes[INDEX2(1, i, NN)],
1631                  elements->Nodes[INDEX2(2, i, NN)],
1632                  elements->Nodes[INDEX2(3, i, NN)],
1633                  elements->Nodes[INDEX2(4, i, NN)],
1634                  elements->Nodes[INDEX2(5, i, NN)],
1635                  elements->Nodes[INDEX2(6, i, NN)],
1636                  elements->Nodes[INDEX2(7, i, NN)],
1637                  elements->Nodes[INDEX2(8, i, NN)],
1638                  elements->Nodes[INDEX2(9, i, NN)],
1639                  elements->Nodes[INDEX2(10, i, NN)],
1640                  elements->Nodes[INDEX2(11, i, NN)],
1641                  elements->Nodes[INDEX2(16, i, NN)],
1642                  elements->Nodes[INDEX2(17, i, NN)],
1643                  elements->Nodes[INDEX2(18, i, NN)],
1644                  elements->Nodes[INDEX2(19, i, NN)],
1645                  elements->Nodes[INDEX2(12, i, NN)],
1646                  elements->Nodes[INDEX2(13, i, NN)],
1647                  elements->Nodes[INDEX2(14, i, NN)],
1648                  elements->Nodes[INDEX2(15, i, NN)]);
1649        }
1650      }
1651      else if (numVTKNodesPerElement!=NN)
1652      {
1653        for (i = 0; i < numCells; i++)
1654        {
1655          for (j = 0; j < numVTKNodesPerElement; j++) fprintf(fileHandle_p,"%d ", elements->Nodes[INDEX2(elements->ReferenceElement->Type->geoNodes[j], i, NN)]);
1656          fprintf(fileHandle_p, "\n");
1657        }
1658      }
1659      else
1660      {
1661        for (i = 0; i < numCells; i++)
1662        {
1663          for (j = 0; j < numVTKNodesPerElement; j++) fprintf(fileHandle_p,"%d ", elements->Nodes[INDEX2(j, i, NN)]);
1664          fprintf(fileHandle_p, "\n");
1665        }
1666      }
1667      fprintf(fileHandle_p, "</DataArray>\n");
1668    
1669      /* write out the DataArray element for the offsets */
1670      fprintf(fileHandle_p, "<DataArray Name=\"offsets\" type=\"Int32\" format=\"ascii\">\n");
1671      for (i=numVTKNodesPerElement; i<=numCells*numVTKNodesPerElement; i+=numVTKNodesPerElement) fprintf(fileHandle_p, "%d\n", i);
1672      fprintf(fileHandle_p, "</DataArray>\n");
1673    
1674      /* write out the DataArray element for the types */
1675      fprintf(fileHandle_p, "<DataArray Name=\"types\" type=\"UInt8\" format=\"ascii\">\n");
1676      for (i=0; i<numCells; i++) fprintf(fileHandle_p, "%d\n", cellType);
1677      fprintf(fileHandle_p, "</DataArray>\n");
1678    
1679      /* finish off the <Cells> element */
1680      fprintf(fileHandle_p, "</Cells>\n");
1681    
1682      /* cell data */
1683      if (write_celldata)
1684      {
1685        /* mark the active data arrays */
1686        bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
1687        fprintf(fileHandle_p, "<CellData");
1688        for (i_data =0 ;i_data<num_data;++i_data)
1689        {
1690          if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data])
1691          {
1692            /* if the rank == 0:   --> scalar data
1693            * if the rank == 1:   --> vector data
1694            * if the rank == 2:   --> tensor data
1695            */
1696            switch(getDataPointRank(data_pp[i_data]))
1697            {
1698            case 0:
1699              if (! set_scalar)
1700              {
1701                fprintf(fileHandle_p," Scalars=\"%s\"",names_p[i_data]);
1702                set_scalar=TRUE;
1703              }
1704              break;
1705            case 1:
1706              if (! set_vector)
1707              {
1708                fprintf(fileHandle_p," Vectors=\"%s\"",names_p[i_data]);
1709                set_vector=TRUE;
1710              }
1711              break;
1712            case 2:
1713              if (! set_tensor)
1714              {
1715                fprintf(fileHandle_p," Tensors=\"%s\"",names_p[i_data]);
1716                set_tensor=TRUE;
1717              }
1718              break;
1719            default:
1720              sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);
1721              Finley_setError(VALUE_ERROR,error_msg);
1722              fclose(fileHandle_p);
1723              return;
1724            }
1725          }
1726        }
1727        fprintf(fileHandle_p, ">\n");
1728        /* write the arrays */
1729        for (i_data =0 ;i_data<num_data;++i_data)
1730        {
1731          if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data])
1732          {
1733            numPointsPerSample = elements->ReferenceElement->numQuadNodes;
1734            rank = getDataPointRank(data_pp[i_data]);
1735            nComp = getDataPointSize(data_pp[i_data]);
1736            nCompReqd=1;   /* the number of components required by vtk */
1737            shape=0;
1738            if (rank == 0)
1739            {
1740              nCompReqd = 1;
1741            }
1742            else if (rank == 1)
1743            {
1744              shape=getDataPointShape(data_pp[i_data], 0);
1745              if  (shape>3)
1746              {
1747                Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");
1748                fclose(fileHandle_p);
1749                return;
1750              }
1751              nCompReqd = 3;
1752            }
1753            else
1754            {
1755              shape=getDataPointShape(data_pp[i_data], 0);
1756              if  (shape>3 || shape != getDataPointShape(data_pp[i_data], 1))
1757              {
1758                Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
1759                fclose(fileHandle_p);
1760                return;
1761              }
1762              nCompReqd = 9;
1763            }
1764            fprintf(fileHandle_p, "<DataArray Name=\"%s\" type=\"Float32\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);
1765    
1766            double sampleAvg[nComp];
1767            for (i=0; i<numCells; i++)
1768            {
1769              values = getSampleData(data_pp[i_data], i);
1770              /* averaging over the number of points in the sample */
1771              for (k=0; k<nComp; k++)
1772              {
1773                rtmp = 0.;
1774                for (j=0; j<numPointsPerSample; j++) rtmp += values[INDEX2(k,j,nComp)];
1775                sampleAvg[k] = rtmp/numPointsPerSample;
1776              }
1777              /* if the number of required components is more than the number
1778              * of actual components, pad with zeros
1779              */
1780              /* probably only need to get shape of first element */
1781              /* write the data different ways for scalar, vector and tensor */
1782              if (nCompReqd == 1)
1783              {
1784                fprintf(fileHandle_p, " %e", sampleAvg[0]);
1785              }
1786              else if (nCompReqd == 3)
1787              {
1788                /* write out the data */
1789                for (m=0; m<shape; m++) fprintf(fileHandle_p, " %e", sampleAvg[m]);
1790                for (m=0; m<nCompReqd-shape; m++) fprintf(fileHandle_p, " %e", 0.);
1791              }
1792              else if (nCompReqd == 9)
1793              {
1794                /* tensor data, so have a 3x3 matrix to output as a row
1795                * of 9 data points */
1796                count = 0;
1797                for (m=0; m<shape; m++)
1798                {
1799                  for (n=0; n<shape; n++)
1800                  {
1801                    fprintf(fileHandle_p, " %e", sampleAvg[count]);
1802                    count++;
1803                  }
1804                  for (n=0; n<3-shape; n++) fprintf(fileHandle_p, " %e", 0.);
1805                }
1806                for (m=0; m<3-shape; m++)
1807                  for (n=0; n<3; n++) fprintf(fileHandle_p, " %e", 0.);
1808              }
1809              fprintf(fileHandle_p, "\n");
1810            }
1811            fprintf(fileHandle_p, "</DataArray>\n");
1812          }
1813        }
1814        fprintf(fileHandle_p, "</CellData>\n");
1815      }
1816      /* point data */
1817      if (write_pointdata)
1818      {
1819        /* mark the active data arrays */
1820        bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
1821        fprintf(fileHandle_p, "<PointData");
1822        for (i_data =0 ;i_data<num_data;++i_data)
1823        {
1824          if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])
1825          {
1826            /* if the rank == 0:   --> scalar data
1827            * if the rank == 1:   --> vector data
1828            * if the rank == 2:   --> tensor data
1829            */
1830            switch(getDataPointRank(data_pp[i_data]))
1831            {
1832            case 0:
1833              if (! set_scalar)
1834              {
1835                fprintf(fileHandle_p," Scalars=\"%s\"",names_p[i_data]);
1836                set_scalar=TRUE;
1837              }
1838              break;
1839            case 1:
1840              if (! set_vector)
1841              {
1842                fprintf(fileHandle_p," Vectors=\"%s\"",names_p[i_data]);
1843                set_vector=TRUE;
1844              }
1845              break;
1846            case 2:
1847              if (! set_tensor)
1848              {
1849                fprintf(fileHandle_p," Tensors=\"%s\"",names_p[i_data]);
1850                set_tensor=TRUE;
1851              }
1852              break;
1853            default:
1854              sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);
1855              Finley_setError(VALUE_ERROR,error_msg);
1856              fclose(fileHandle_p);
1857              return;
1858            }
1859          }
1860        }
1861        fprintf(fileHandle_p, ">\n");
1862        /* write the arrays */
1863        for (i_data =0 ;i_data<num_data;++i_data)
1864        {
1865          if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])
1866          {
1867            numPointsPerSample = elements->ReferenceElement->numQuadNodes;
1868            rank = getDataPointRank(data_pp[i_data]);
1869            nComp = getDataPointSize(data_pp[i_data]);
1870            nCompReqd=1;   /* the number of components required by vtk */
1871            shape=0;
1872            if (rank == 0)
1873            {
1874              nCompReqd = 1;
1875            }
1876            else if (rank == 1)
1877            {
1878              shape=getDataPointShape(data_pp[i_data], 0);
1879              if  (shape>3)
1880              {
1881                Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");
1882                fclose(fileHandle_p);
1883                return;
1884              }
1885              nCompReqd = 3;
1886            }
1887            else
1888            {
1889              shape=getDataPointShape(data_pp[i_data], 0);
1890              if  (shape>3 || shape != getDataPointShape(data_pp[i_data], 1))
1891              {
1892                Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
1893                fclose(fileHandle_p);
1894                return;
1895              }
1896              nCompReqd = 9;
1897            }
1898            fprintf(fileHandle_p, "<DataArray Name=\"%s\" type=\"Float32\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);
1899            /* write out the data */
1900            /* if the number of required components is more than the number
1901            * of actual components, pad with zeros
1902            */
1903            bool_t do_write=TRUE;
1904            for (i=0; i<mesh_p->Nodes->numNodes; i++)
1905            {
1906              if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1907              {
1908                if (mesh_p->Nodes->toReduced[i]>=0)
1909                {
1910                  switch(getFunctionSpaceType(data_pp[i_data]))
1911                  {
1912                  case FINLEY_DEGREES_OF_FREEDOM:
1913                    values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);
1914                    break;
1915                  case FINLEY_REDUCED_DEGREES_OF_FREEDOM:
1916                    values = getSampleData(data_pp[i_data],mesh_p->Nodes->reducedDegreeOfFreedom[i]);
1917                    break;
1918                  case FINLEY_NODES:
1919                    values = getSampleData(data_pp[i_data],i);
1920                    break;
1921                  }
1922                  do_write=TRUE;
1923                }
1924                else
1925                {
1926                  do_write=FALSE;
1927                }
1928              }
1929              else
1930              {
1931                do_write=TRUE;
1932                switch(getFunctionSpaceType(data_pp[i_data]))
1933                {
1934                case FINLEY_DEGREES_OF_FREEDOM:
1935                  values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);
1936                  break;
1937                case FINLEY_NODES:
1938                  values = getSampleData(data_pp[i_data],i);
1939                  break;
1940                }
1941              }
1942              /* write the data different ways for scalar, vector and tensor */
1943              if (do_write)
1944              {
1945                if (nCompReqd == 1)
1946                {
1947                  fprintf(fileHandle_p, " %e", values[0]);
1948                }
1949                else if (nCompReqd == 3)
1950                {
1951                  for (m=0; m<shape; m++) fprintf(fileHandle_p, " %e", values[m]);
1952                  for (m=0; m<nCompReqd-shape; m++) fprintf(fileHandle_p, " %e", 0.);
1953                }
1954                else if (nCompReqd == 9)
1955                {
1956                  /* tensor data, so have a 3x3 matrix to output as a row
1957                  * of 9 data points */
1958                  count = 0;
1959                  for (m=0; m<shape; m++)
1960                  {
1961                    for (n=0; n<shape; n++)
1962                    {
1963                      fprintf(fileHandle_p, " %e", values[count]);
1964                      count++;
1965                    }
1966                    for (n=0; n<3-shape; n++) fprintf(fileHandle_p, " %e", 0.);
1967                  }
1968                  for (m=0; m<3-shape; m++)
1969                    for (n=0; n<3; n++) fprintf(fileHandle_p, " %e", 0.);
1970                }
1971                fprintf(fileHandle_p, "\n");
1972              }
1973            }
1974            fprintf(fileHandle_p, "</DataArray>\n");
1975          }
1976        }
1977        fprintf(fileHandle_p, "</PointData>\n");
1978      }
1979    /* finish off the piece */    /* finish off the piece */
1980    fprintf(fileHandle_p, "</Piece>\n");    fprintf(fileHandle_p, "</Piece>\n");
1981    
# Line 614  void Finley_Mesh_saveVTK(const char * fi Line 1986  void Finley_Mesh_saveVTK(const char * fi
1986    fclose(fileHandle_p);    fclose(fileHandle_p);
1987    return;    return;
1988  }  }
1989    #endif
1990    
 /*  
  * $Log$  
  * Revision 1.4  2005/05/06 04:26:15  jgs  
  * Merge of development branch back to main trunk on 2005-05-06  
  *  
  * Revision 1.1.2.6  2005/05/06 01:17:19  cochrane  
  * Fixed incorrect reporting of number of components in PointData arrays for  
  * vector data.  
  *  
  * Revision 1.1.2.5  2005/05/05 05:38:44  cochrane  
  * Improved formatting of VTK file output.  
  *  
  * Revision 1.1.2.4  2005/02/22 10:03:54  cochrane  
  * Implementation of writing of vtk xml files from finley.  This function will  
  * require more testing, but on the cases that I have tried (and with the help  
  * of Lutz and mayavi), it looks like it's producing the correct output.  Testing  
  * with more realistic data would be good.  I'm at least confident enough to  
  * commit my changes.  
  *  
  * Revision 1.1.2.3  2005/02/17 05:53:26  gross  
  * some bug in saveDX fixed: in fact the bug was in  
  * DataC/getDataPointShape  
  *  
  * Revision 1.1.2.2  2005/02/10 01:34:22  cochrane  
  * Quick fix to make sure that saveVTK compiles so that finley is still buildable.  Apologies to those this has affected.  
  *  
  * Revision 1.1.2.1  2005/02/09 06:53:15  cochrane  
  * Initial import to repository.  This is the file to implement saving finley/escript meshes out to vtk formatted files.  It is basically just a hack of the opendx equivalent, with a lot of the opendx stuff still in the file, so it doesn't actually work just yet, but it probably needs to be added to the cvs.  
  *  
  * Revision 1.1.1.1  2004/10/26 06:53:57  jgs  
  * initial import of project esys2  
  *  
  * Revision 1.1  2004/07/27 08:27:11  gross  
  * Finley: saveDX added: now it is possible to write data on boundary and contact elements  
  *  
  */  

Legend:
Removed from v.121  
changed lines
  Added in v.796

  ViewVC Help
Powered by ViewVC 1.1.26