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

  ViewVC Help
Powered by ViewVC 1.1.26