/[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 110 by jgs, Mon Feb 14 04:14:42 2005 UTC trunk/finley/src/Mesh_saveVTK.c revision 818 by dhawcroft, Sun Aug 27 11:10:34 2006 UTC
# Line 1  Line 1 
1  /* $Id$ */  /*
2    ************************************************************
3    *          Copyright 2006 by ACcESS MNRF                   *
4    *                                                          *
5    *              http://www.access.edu.au                    *
6    *       Primary Business: Queensland, Australia            *
7    *  Licensed under the Open Software License version 3.0    *
8    *     http://www.opensource.org/licenses/osl-3.0.php       *
9    *                                                          *
10    ************************************************************
11    */
12    
13    
14  /**************************************************************/  /**************************************************************/
15    
# Line 6  Line 17 
17    
18  /**************************************************************/  /**************************************************************/
19    
 /*   Copyrights by ACcESS, Australia 2004 */  
20  /*   Author: Paul Cochrane, cochrane@esscc.uq.edu.au */  /*   Author: Paul Cochrane, cochrane@esscc.uq.edu.au */
21    /*   MPI-IO version: Derick Hawcroft, d.hawcroft@uq.edu.au         */
22    
23    /*   Version: $Id$ */
24    
25  /**************************************************************/  /**************************************************************/
26    
27  #include "Finley.h"  
 #include "Common.h"  
28  #include "Mesh.h"  #include "Mesh.h"
29  #include "escript/Data/DataC.h"  #include "vtkCellType.h"  /* copied from vtk source directory !!! */
30    
31    /*
32     MPI version notes:
33    
34     ******************************************************************************
35     ***                                                                       ****
36     *** 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;
81    
82      // Local element info (for debugging)
83      size_t numLocalCells,
84      numInternalCells,
85      numBoundaryCells;
86    
87      int rank;
88    
89      int amode = MPI_MODE_CREATE | MPI_MODE_WRONLY |  MPI_MODE_SEQUENTIAL;
90    
91      comm   = mesh_p->Nodes->MPIInfo->comm;
92      myRank = mesh_p->Nodes->MPIInfo->rank;
93      gsize  = mesh_p->Nodes->MPIInfo->size;
94    
95      MPI_File fh;
96      MPI_Status status;
97      MPI_Request req;
98      MPI_Info infoHints;
99    
100      char error_msg[LenErrorMsg_MAX];
101    
102      int i_data;
103    
104      int nodetype=FINLEY_DEGREES_OF_FREEDOM;
105      int elementtype=FINLEY_UNKNOWN;
106      bool_t isCellCentered[num_data],write_celldata=FALSE,write_pointdata=FALSE;
107    
108      ElementTypeId TypeId;
109    
110      int numVTKNodesPerElement;
111      int cellType;
112      char elemTypeStr[32];
113    
114      Finley_ElementFile* elements=NULL;
115    
116    
117      // Local node info
118      int numInternalNodes,
119      numLocalNodes,
120      numBoundaryNodes,
121      localDOF;  // equals to  (DOF of Internal Nodes) +  (DOF of Boundary Nodes) of local domain
122    
123    
124      nDim  = mesh_p->Nodes->numDim;
125    
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      Finley_NodeDistribution *dist;
291      if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
292      {
293        dist = mesh_p->Nodes->reducedDegreeOfFreedomDistribution;
294      }
295      else
296      {
297        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;
443      }
444    
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=\"Float32\" 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=\"Float32\" 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                rtmp = 0.;
870                for (j=0; j<numPointsPerSample; j++) rtmp += values[INDEX2(n,j,nComp)];
871                sampleAvg[k] = rtmp/numPointsPerSample;
872              }
873              // if the number of required components is more than the number
874              // of actual components, pad with zeros
875    
876              // probably only need to get shape of first element
877              // write the data different ways for scalar, vector and tensor
878              if (nCompReqd == 1)
879              {
880                sprintf(tmpbuf, " %e", sampleAvg[0]);
881                __STRCAT(largebuf,tmpbuf,tsz)
882              }
883              else if (nCompReqd == 3)
884              {
885                // write out the data
886                for (m=0; m<shape; m++)
887                {
888                  sprintf(tmpbuf, " %e", sampleAvg[m]);
889                  __STRCAT(largebuf,tmpbuf,tsz)
890                }
891                for (m=0; m<nCompReqd-shape; m++)
892                {
893                  __STRCAT(largebuf,zero,tsz)
894                }
895              }
896              else if (nCompReqd == 9)
897              {
898                // tensor data, so have a 3x3 matrix to output as a row
899                // of 9 data points
900                count = 0;
901                for (m=0; m<shape; m++)
902                {
903                  for (n=0; n<shape; n++)
904                  {
905                    sprintf(tmpbuf, " %e", sampleAvg[count]);
906                    __STRCAT(largebuf,tmpbuf,tsz)
907                    count++;
908                  }
909                  for (n=0; n<3-shape; n++)
910                  {
911                    __STRCAT(largebuf,zero,tsz)
912                  }
913                }
914                for (m=0; m<3-shape; m++)
915                  for (n=0; n<3; n++)
916                  {
917                    __STRCAT(largebuf,zero,tsz)
918                  }
919              }
920              __STRCAT(largebuf,newline,tsz)
921            }
922            largebuf[tsz] = '\0';
923            MPI_File_write_ordered(fh,largebuf,tsz,MPI_CHAR,&status);
924            MEMFREE(largebuf);
925            if( myRank == 0)
926            {
927              char *tag = "</DataArray>\n";
928              MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);
929              MPI_Wait(&req,&status);
930            }
931    
932          }
933        }
934        // closing celldata tag
935        if(myRank == 0)
936        {
937          char* tag =  "</CellData>\n";
938          MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);
939          MPI_Wait(&req,&status);
940        }
941    
942        MPIO_DEBUG(" Done Writing Cell Data ")
943      }
944    
945    
946      // Write Point Data Header Tags
947      if( myRank == 0)
948      {
949        char header[600];
950        char tmpBuf[50];
951    
952        if (write_pointdata)
953        {
954          MPIO_DEBUG(" Writing Pointdata... ")
955          // mark the active data arrays
956          bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
957          sprintf(header, "<PointData");
958          for (i_data =0 ;i_data<num_data;++i_data)
959          {
960            if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])
961            {
962              // if the rank == 0:   --> scalar data
963              // if the rank == 1:   --> vector data
964              // if the rank == 2:   --> tensor data
965    
966              switch(getDataPointRank(data_pp[i_data]))
967              {
968              case 0:
969                if (! set_scalar)
970                {
971                  sprintf(tmpBuf," Scalars=\"%s\"",names_p[i_data]);
972                  strcat(header,tmpBuf);
973                  set_scalar=TRUE;
974                }
975                break;
976              case 1:
977                if (! set_vector)
978                {
979                  sprintf(tmpBuf," Vectors=\"%s\"",names_p[i_data]);
980                  strcat(header,tmpBuf);
981                  set_vector=TRUE;
982                }
983                break;
984              case 2:
985                if (! set_tensor)
986                {
987                  sprintf(tmpBuf," Tensors=\"%s\"",names_p[i_data]);
988                  strcat(header,tmpBuf);
989                  set_tensor=TRUE;
990                }
991                break;
992              default:
993                sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);
994                Finley_setError(VALUE_ERROR,error_msg);
995                return;
996              }
997            }
998          }
999          strcat(header, ">\n");
1000          MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);
1001          MPI_Wait(&req,&status);
1002        }
1003      }
1004    
1005      // write actual data
1006      if(write_pointdata)
1007      {
1008        for (i_data =0 ;i_data<num_data;++i_data)
1009        {
1010          if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])
1011          {
1012            numPointsPerSample = elements->ReferenceElement->numQuadNodes;
1013            rank = getDataPointRank(data_pp[i_data]);
1014            nComp = getDataPointSize(data_pp[i_data]);
1015            nCompReqd=1;   // the number of components required by vtk
1016            shape=0;
1017            if (rank == 0)
1018            {
1019              nCompReqd = 1;
1020            }
1021            else if (rank == 1)
1022            {
1023              shape=getDataPointShape(data_pp[i_data], 0);
1024              if  (shape>3)
1025              {
1026                Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");
1027                return;
1028              }
1029              nCompReqd = 3;
1030            }
1031            else
1032            {
1033              shape=getDataPointShape(data_pp[i_data], 0);
1034              if  (shape>3 || shape != getDataPointShape(data_pp[i_data], 1))
1035              {
1036                Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
1037                return;
1038              }
1039              nCompReqd = 9;
1040            }
1041    
1042            if( myRank == 0)
1043            {
1044              char header[250];
1045              sprintf(header, "<DataArray Name=\"%s\" type=\"Float32\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);
1046              MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);
1047              MPI_Wait(&req,&status);
1048            }
1049            // write out the data
1050            // if the number of required components is more than the number
1051            // of actual components, pad with zeros
1052    
1053            char tmpbuf[15];
1054            char* largebuf = MEMALLOC(nCompReqd*numLocalNodes*15 + numLocal*nCompReqd + nCompReqd + 15,char);
1055            largebuf[0] = '\0';
1056            bool_t do_write=TRUE;
1057            size_t tsz = 0;
1058    
1059            for(k=0;k < nodeCache.size;k++)
1060            {
1061              i = nodeCache.values[k];
1062    
1063              if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1064              {
1065                if (mesh_p->Nodes->toReduced[i]>=0)
1066                {
1067                  switch(getFunctionSpaceType(data_pp[i_data]))
1068                  {
1069                  case FINLEY_DEGREES_OF_FREEDOM:
1070                    values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);
1071                    break;
1072                  case FINLEY_REDUCED_DEGREES_OF_FREEDOM:
1073                    values = getSampleData(data_pp[i_data],mesh_p->Nodes->reducedDegreeOfFreedom[i]);
1074                    break;
1075                  case FINLEY_NODES:
1076                    values = getSampleData(data_pp[i_data],i);
1077                    break;
1078                  }
1079                  do_write=TRUE;
1080                }
1081                else
1082                {
1083                  do_write=FALSE;
1084                }
1085              }
1086              else
1087              {
1088                do_write=TRUE;
1089                switch(getFunctionSpaceType(data_pp[i_data]))
1090                {
1091                case FINLEY_DEGREES_OF_FREEDOM:
1092                  values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);
1093                  break;
1094                case FINLEY_NODES:
1095                  values = getSampleData(data_pp[i_data],i);
1096                  break;
1097                }
1098              }
1099              // write the data different ways for scalar, vector and tensor
1100              if (do_write)
1101              {
1102                if (nCompReqd == 1)
1103                {
1104                  sprintf(tmpbuf," %e", values[0]);
1105                  __STRCAT(largebuf,tmpbuf,tsz)
1106                }
1107                else if (nCompReqd == 3)
1108                {
1109                  for (m=0; m<shape; m++)
1110                  {
1111    
1112                    sprintf(tmpbuf," %e",values[m]);
1113                    __STRCAT(largebuf,tmpbuf,tsz)
1114                  }
1115                  for (m=0; m<nCompReqd-shape; m++)
1116                  {
1117                    __STRCAT(largebuf,zero,tsz)
1118                  }
1119                }
1120                else if (nCompReqd == 9)
1121                {
1122                  // tensor data, so have a 3x3 matrix to output as a row
1123                  //  of 9 data points
1124                  count = 0;
1125                  for (m=0; m<shape; m++)
1126                  {
1127                    for (n=0; n<shape; n++)
1128                    {
1129                      sprintf(tmpbuf," %e", values[count]);
1130                      __STRCAT(largebuf,tmpbuf,tsz)
1131                      count++;
1132                    }
1133                    for (n=0; n<3-shape; n++)
1134                    {
1135                      __STRCAT(largebuf,zero,tsz)
1136                    }
1137                  }
1138                  for (m=0; m<3-shape; m++)
1139                  {
1140                    for (n=0; n<3; n++)
1141                    {
1142                      __STRCAT(largebuf,zero,tsz)
1143                    }
1144                  }
1145                }
1146                __STRCAT(largebuf,newline,tsz)
1147              }
1148    
1149  void Finley_Mesh_saveVTK(const char * filename_p, Finley_Mesh *mesh_p, escriptDataC* data_p) {          }
1150            // Write out local data
1151    
1152            largebuf[tsz] = '\0';
1153            MPI_File_write_ordered(fh,largebuf,tsz,MPI_CHAR,&status);
1154            MEMFREE(largebuf);
1155            if( myRank == 0)
1156            {
1157              char *tag = "</DataArray>\n";
1158              MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);
1159              MPI_Wait(&req,&status);
1160            }
1161          }
1162        }
1163        // Finish off with closing tag
1164        if(myRank == 0)
1165        {
1166          char* tag =  "</PointData>\n";
1167          MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);
1168          MPI_Wait(&req,&status);
1169        }
1170        MPIO_DEBUG(" Done Writing Pointdata ")
1171      }
1172      // end write_pointdata
1173    
1174      // tag and bag...  
1175      if (myRank == 0)
1176      {
1177        char *footer = "</Piece>\n</UnstructuredGrid>\n</VTKFile>";
1178        MPI_File_iwrite_shared(fh,footer,strlen(footer),MPI_CHAR,&req);
1179        MPI_Wait(&req,&status);
1180      }
1181    
1182      MEMFREE(nodesGlobal);
1183      MEMFREE(nodeCache.values);
1184      MEMFREE(elementCache.values);
1185    #ifdef MPIO_HINTS
1186      MPI_Info_free(&infoHints);
1187    #undef MPIO_HINTS
1188    #endif
1189      MPI_File_close(&fh);
1190      MPIO_DEBUG(" ***** Exit saveVTK ***** ")
1191    #undef __STRCAT
1192    }
1193    
1194    #undef MPIO_DEBUG
1195    #else
1196    
1197    
1198    
1199    
1200    void Finley_Mesh_saveVTK(const char * filename_p, Finley_Mesh *mesh_p, const dim_t num_data,char* *names_p, escriptDataC* *data_pp)
1201    {
1202      char error_msg[LenErrorMsg_MAX];
1203    /* if there is no mesh we just return */    /* if there is no mesh we just return */
1204    if (mesh_p==NULL) return;    if (mesh_p==NULL) return;
1205    /* some tables needed for reordering */  
1206    int resort[6][9]={    int i, j, k, numVTKNodesPerElement,i_data,m, count, n, rank,shape, numPoints, cellType, numCells,
1207                      {0,1},             /* line */    nDim, numPointsPerSample, nComp, nCompReqd;
1208                      {0,1,2},           /* triangle */  
1209                      {0,1,2,3},         /* tetrahedron */    index_t j2;
1210                      {0,3,1,2},         /* quadrilateral */    double* values, rtmp;
                     {3,0,7,4,2,1,6,5}, /* hexahedron */  
                    };  
   Finley_ElementFile* elements=NULL;  
1211    char elemTypeStr[32];    char elemTypeStr[32];
1212    int i,j,k,numVTKNodesPerElement,*resortIndex,isCellCentered,nodetype;  
1213    double* values,rtmp;    /* open the file and check handle */
1214    int nDim = mesh_p->Nodes->numDim;  
   /* open the file  and check handle */  
1215    FILE * fileHandle_p = fopen(filename_p, "w");    FILE * fileHandle_p = fopen(filename_p, "w");
1216    if (fileHandle_p==NULL) {    if (fileHandle_p==NULL)
1217           Finley_ErrorCode=IO_ERROR;    {
1218           sprintf(Finley_ErrorMsg,"File %s could not be opened for writing.",filename_p);      sprintf(error_msg, "saveVTK: File %s could not be opened for writing.", filename_p);
1219           return;      Finley_setError(IO_ERROR,error_msg);
1220        return;
1221      }
1222      /* find the mesh type to be written */
1223      int nodetype=FINLEY_DEGREES_OF_FREEDOM;
1224      int elementtype=FINLEY_UNKNOWN;
1225      bool_t isCellCentered[num_data],write_celldata=FALSE,write_pointdata=FALSE;
1226      for (i_data=0;i_data<num_data;++i_data)
1227      {
1228        if (! isEmpty(data_pp[i_data]))
1229        {
1230          switch(getFunctionSpaceType(data_pp[i_data]))
1231          {
1232          case FINLEY_DEGREES_OF_FREEDOM:
1233            nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1234            if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)
1235            {
1236              elementtype=FINLEY_ELEMENTS;
1237            }
1238            else
1239            {
1240              Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1241              fclose(fileHandle_p);
1242              return;
1243            }
1244            isCellCentered[i_data]=FALSE;
1245            break;
1246          case FINLEY_REDUCED_DEGREES_OF_FREEDOM:
1247            nodetype = FINLEY_REDUCED_DEGREES_OF_FREEDOM;
1248            if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)
1249            {
1250              elementtype=FINLEY_ELEMENTS;
1251            }
1252            else
1253            {
1254              Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1255              fclose(fileHandle_p);
1256              return;
1257            }
1258            isCellCentered[i_data]=FALSE;
1259            break;
1260          case FINLEY_NODES:
1261            nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1262            if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)
1263            {
1264              elementtype=FINLEY_ELEMENTS;
1265            }
1266            else
1267            {
1268              Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1269              fclose(fileHandle_p);
1270              return;
1271            }
1272            isCellCentered[i_data]=FALSE;
1273            break;
1274          case FINLEY_ELEMENTS:
1275            nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1276            if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)
1277            {
1278              elementtype=FINLEY_ELEMENTS;
1279            }
1280            else
1281            {
1282              Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1283              fclose(fileHandle_p);
1284              return;
1285            }
1286            isCellCentered[i_data]=TRUE;
1287            break;
1288          case FINLEY_FACE_ELEMENTS:
1289            nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1290            if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_FACE_ELEMENTS)
1291            {
1292              elementtype=FINLEY_FACE_ELEMENTS;
1293            }
1294            else
1295            {
1296              Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1297              fclose(fileHandle_p);
1298              return;
1299            }
1300            isCellCentered[i_data]=TRUE;
1301            break;
1302          case FINLEY_POINTS:
1303            nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1304            if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_POINTS)
1305            {
1306              elementtype=FINLEY_POINTS;
1307            }
1308            else
1309            {
1310              Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1311              fclose(fileHandle_p);
1312              return;
1313            }
1314            isCellCentered[i_data]=TRUE;
1315            break;
1316          case FINLEY_CONTACT_ELEMENTS_1:
1317            nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1318            if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1)
1319            {
1320              elementtype=FINLEY_CONTACT_ELEMENTS_1;
1321            }
1322            else
1323            {
1324              Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1325              fclose(fileHandle_p);
1326              return;
1327            }
1328            isCellCentered[i_data]=TRUE;
1329            break;
1330          case FINLEY_CONTACT_ELEMENTS_2:
1331            nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1332            if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1)
1333            {
1334              elementtype=FINLEY_CONTACT_ELEMENTS_1;
1335            }
1336            else
1337            {
1338              Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1339              fclose(fileHandle_p);
1340              return;
1341            }
1342            isCellCentered[i_data]=TRUE;
1343            break;
1344          default:
1345            sprintf(error_msg,"saveVTK: Finley does not know anything about function space type %d",getFunctionSpaceType(data_pp[i_data]));
1346            Finley_setError(TYPE_ERROR,error_msg);
1347            fclose(fileHandle_p);
1348            return;
1349          }
1350          if (isCellCentered[i_data])
1351          {
1352            write_celldata=TRUE;
1353          }
1354          else
1355          {
1356            write_pointdata=TRUE;
1357          }
1358        }
1359      }
1360      /* select nomber of points and the mesh component */
1361      numPoints = mesh_p->Nodes->numNodes;
1362      if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1363      {
1364        numPoints = mesh_p->Nodes->reducedNumNodes;
1365      }
1366      else
1367      {
1368        numPoints = mesh_p->Nodes->numNodes;
1369      }
1370      if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;
1371      Finley_ElementFile* elements=NULL;
1372      switch(elementtype)
1373      {
1374      case FINLEY_ELEMENTS:
1375        elements=mesh_p->Elements;
1376        break;
1377      case FINLEY_FACE_ELEMENTS:
1378        elements=mesh_p->FaceElements;
1379        break;
1380      case FINLEY_POINTS:
1381        elements=mesh_p->Points;
1382        break;
1383      case FINLEY_CONTACT_ELEMENTS_1:
1384        elements=mesh_p->ContactElements;
1385        break;
1386      }
1387      if (elements==NULL)
1388      {
1389        Finley_setError(SYSTEM_ERROR,"saveVTK: undefined element file");
1390        fclose(fileHandle_p);
1391        return;
1392      }
1393      /* map finley element type to VTK element type */
1394      numCells = elements->numElements;
1395      ElementTypeId TypeId;
1396      if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1397      {
1398        TypeId = elements->LinearReferenceElement->Type->TypeId;
1399      }
1400      else
1401      {
1402        TypeId = elements->ReferenceElement->Type->TypeId;
1403      }
1404    
1405      switch(TypeId)
1406      {
1407      case Point1:
1408      case Line2Face:
1409      case Line3Face:
1410      case Point1_Contact:
1411      case Line2Face_Contact:
1412      case Line3Face_Contact:
1413        cellType = VTK_VERTEX;
1414        numVTKNodesPerElement = 1;
1415        strcpy(elemTypeStr, "VTK_VERTEX");
1416        break;
1417    
1418      case Line2:
1419      case Tri3Face:
1420      case Rec4Face:
1421      case Line2_Contact:
1422      case Tri3_Contact:
1423      case Tri3Face_Contact:
1424      case Rec4Face_Contact:
1425        cellType = VTK_LINE;
1426        numVTKNodesPerElement = 2;
1427        strcpy(elemTypeStr, "VTK_LINE");
1428        break;
1429    
1430      case Tri3:
1431      case Tet4Face:
1432      case Tet4Face_Contact:
1433        cellType = VTK_TRIANGLE;
1434        numVTKNodesPerElement = 3;
1435        strcpy(elemTypeStr, "VTK_TRIANGLE");
1436        break;
1437    
1438      case Rec4:
1439      case Hex8Face:
1440      case Rec4_Contact:
1441      case Hex8Face_Contact:
1442        cellType = VTK_QUAD;
1443        numVTKNodesPerElement = 4;
1444        strcpy(elemTypeStr, "VTK_QUAD");
1445        break;
1446    
1447      case Tet4:
1448        cellType = VTK_TETRA;
1449        numVTKNodesPerElement = 4;
1450        strcpy(elemTypeStr, "VTK_TETRA");
1451        break;
1452    
1453      case Hex8:
1454        cellType = VTK_HEXAHEDRON;
1455        numVTKNodesPerElement = 8;
1456        strcpy(elemTypeStr, "VTK_HEXAHEDRON");
1457        break;
1458    
1459      case Line3:
1460      case Tri6Face:
1461      case Rec8Face:
1462      case Line3_Contact:
1463      case Tri6Face_Contact:
1464      case Rec8Face_Contact:
1465        cellType = VTK_QUADRATIC_EDGE;
1466        numVTKNodesPerElement = 3;
1467        strcpy(elemTypeStr, "VTK_QUADRATIC_EDGE");
1468        break;
1469    
1470      case Tri6:
1471      case Tet10Face:
1472      case Tri6_Contact:
1473      case Tet10Face_Contact:
1474        cellType = VTK_QUADRATIC_TRIANGLE;
1475        numVTKNodesPerElement = 6;
1476        strcpy(elemTypeStr, "VTK_QUADRATIC_TRIANGLE");
1477        break;
1478    
1479      case Rec8:
1480      case Hex20Face:
1481      case Rec8_Contact:
1482      case Hex20Face_Contact:
1483        cellType = VTK_QUADRATIC_QUAD;
1484        numVTKNodesPerElement = 8;
1485        strcpy(elemTypeStr, "VTK_QUADRATIC_QUAD");
1486        break;
1487    
1488      case Tet10:
1489        cellType = VTK_QUADRATIC_TETRA;
1490        numVTKNodesPerElement = 10;
1491        strcpy(elemTypeStr, "VTK_QUADRATIC_TETRA");
1492        break;
1493    
1494      case Hex20:
1495        cellType = VTK_QUADRATIC_HEXAHEDRON;
1496        numVTKNodesPerElement = 20;
1497        strcpy(elemTypeStr, "VTK_QUADRATIC_HEXAHEDRON");
1498        break;
1499    
1500      default:
1501        sprintf(error_msg, "saveVTK: Element type %s is not supported by VTK",elements->ReferenceElement->Type->Name);
1502        Finley_setError(VALUE_ERROR,error_msg);
1503        fclose(fileHandle_p);
1504        return;
1505    }    }
1506    /* xml header */    /* xml header */
1507    fprintf(fileHandle_p, "<?xml version=\"1.0\"?>\n");    fprintf(fileHandle_p, "<?xml version=\"1.0\"?>\n");
1508    fprintf(fileHandle_p,    fprintf(fileHandle_p, "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\">\n");
       "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\">\n");  
1509    
1510    /* finley uses an unstructured mesh, so UnstructuredGrid *should* work */    /* finley uses an unstructured mesh, so UnstructuredGrid *should* work */
1511    fprintf(fileHandle_p, "<UnstructuredGrid>\n");    fprintf(fileHandle_p, "<UnstructuredGrid>\n");
1512    
1513    /* is there only one "piece" to the data?? */    /* is there only one "piece" to the data?? */
1514    /* fprintf(fileHandle_p,    fprintf(fileHandle_p, "<Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n",numPoints, numCells);
1515        "<Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n",);    /* now for the points; equivalent to positions section in saveDX() */
1516        */    /* "The points element explicitly defines coordinates for each point
1517      * individually.  It contains one DataArray element describing an array
1518      * with three components per value, each specifying the coordinates of one
1519      * point" - from Vtk User's Guide
1520      */
1521      fprintf(fileHandle_p, "<Points>\n");
1522      /*
1523      * the reason for this if statement is explained in the long comment below
1524      */
1525      nDim = mesh_p->Nodes->numDim;
1526      fprintf(fileHandle_p, "<DataArray NumberOfComponents=\"%d\" type=\"Float32\" format=\"ascii\">\n",MAX(3,nDim));
1527      /* vtk/mayavi doesn't like 2D data, it likes 3D data with a degenerate
1528      * third dimension to handle 2D data (like a sheet of paper).  So, if
1529      * nDim is 2, we have to append zeros to the array to get this third
1530      * dimension, and keep the visualisers happy.
1531      * Indeed, if nDim is less than 3, must pad all empty dimensions, so
1532      * that the total number of dims is 3.
1533      */
1534      if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1535      {
1536        for (i = 0; i < mesh_p->Nodes->numNodes; i++)
1537        {
1538          if (mesh_p->Nodes->toReduced[i]>=0)
1539          {
1540            for (j = 0; j < nDim; j++) fprintf(fileHandle_p, " %e",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);
1541            for (k=0; k<3-nDim; k++) fprintf(fileHandle_p, " %e",0.);
1542            fprintf(fileHandle_p, "\n");
1543          }
1544        }
1545      }
1546      else
1547      {
1548        for (i = 0; i < mesh_p->Nodes->numNodes; i++)
1549        {
1550    
1551          for (j = 0; j < nDim; j++) fprintf(fileHandle_p, " %e",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);
1552          for (k=0; k<3-nDim; k++) fprintf(fileHandle_p, " %e",0.);
1553          fprintf(fileHandle_p, "\n");
1554        }
1555    
1556    /* now for the point data */    }
1557    /* fprintf(fileHandle_p,    fprintf(fileHandle_p, "</DataArray>\n");
1558        "<PointData Scalars=\"%s\" Vectors=\"%s\">\n",);    fprintf(fileHandle_p, "</Points>\n");
       */  
1559    
1560    fprintf(fileHandle_p, "</PointData>\n");    /* write out the DataArray element for the connectivity */
1561    
1562    /* now for the cell data */    int NN = elements->ReferenceElement->Type->numNodes;
1563    /* fprintf(fileHandle_p,    fprintf(fileHandle_p, "<Cells>\n");
1564        "<CellData Scalars=\"%s\" Vectors=\"%s\">\n",);    fprintf(fileHandle_p, "<DataArray Name=\"connectivity\" type=\"Int32\" format=\"ascii\">\n");
       */  
1565    
1566    fprintf(fileHandle_p, "</CellData>\n");    if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1567      {
1568        for (i = 0; i < numCells; i++)
1569        {
1570          for (j = 0; j < numVTKNodesPerElement; j++)
1571            fprintf(fileHandle_p,"%d ",mesh_p->Nodes->toReduced[elements->Nodes[INDEX2(elements->ReferenceElement->Type->linearNodes[j], i, NN)]]);
1572          fprintf(fileHandle_p, "\n");
1573        }
1574      }
1575      else if (VTK_QUADRATIC_HEXAHEDRON==cellType)
1576      {
1577        for (i = 0; i < numCells; i++)
1578        {
1579          fprintf(fileHandle_p,"%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d\n",
1580                  elements->Nodes[INDEX2(0, i, NN)],
1581                  elements->Nodes[INDEX2(1, i, NN)],
1582                  elements->Nodes[INDEX2(2, i, NN)],
1583                  elements->Nodes[INDEX2(3, i, NN)],
1584                  elements->Nodes[INDEX2(4, i, NN)],
1585                  elements->Nodes[INDEX2(5, i, NN)],
1586                  elements->Nodes[INDEX2(6, i, NN)],
1587                  elements->Nodes[INDEX2(7, i, NN)],
1588                  elements->Nodes[INDEX2(8, i, NN)],
1589                  elements->Nodes[INDEX2(9, i, NN)],
1590                  elements->Nodes[INDEX2(10, i, NN)],
1591                  elements->Nodes[INDEX2(11, i, NN)],
1592                  elements->Nodes[INDEX2(16, i, NN)],
1593                  elements->Nodes[INDEX2(17, i, NN)],
1594                  elements->Nodes[INDEX2(18, i, NN)],
1595                  elements->Nodes[INDEX2(19, i, NN)],
1596                  elements->Nodes[INDEX2(12, i, NN)],
1597                  elements->Nodes[INDEX2(13, i, NN)],
1598                  elements->Nodes[INDEX2(14, i, NN)],
1599                  elements->Nodes[INDEX2(15, i, NN)]);
1600        }
1601      }
1602      else if (numVTKNodesPerElement!=NN)
1603      {
1604        for (i = 0; i < numCells; i++)
1605        {
1606          for (j = 0; j < numVTKNodesPerElement; j++) fprintf(fileHandle_p,"%d ", elements->Nodes[INDEX2(elements->ReferenceElement->Type->geoNodes[j], i, NN)]);
1607          fprintf(fileHandle_p, "\n");
1608        }
1609      }
1610      else
1611      {
1612        for (i = 0; i < numCells; i++)
1613        {
1614          for (j = 0; j < numVTKNodesPerElement; j++) fprintf(fileHandle_p,"%d ", elements->Nodes[INDEX2(j, i, NN)]);
1615          fprintf(fileHandle_p, "\n");
1616        }
1617      }
1618      fprintf(fileHandle_p, "</DataArray>\n");
1619    
1620    /* now for the points */    /* write out the DataArray element for the offsets */
1621    fprintf(fileHandle_p, "<Points>\n");    fprintf(fileHandle_p, "<DataArray Name=\"offsets\" type=\"Int32\" format=\"ascii\">\n");
1622    fprintf(fileHandle_p, "<DataArray NumberOfComponents=\"%d\">\n");    for (i=numVTKNodesPerElement; i<=numCells*numVTKNodesPerElement; i+=numVTKNodesPerElement) fprintf(fileHandle_p, "%d\n", i);
1623    fprintf(fileHandle_p, "</DataArray>\n");    fprintf(fileHandle_p, "</DataArray>\n");
   fprintf(fileHandle_p, "</Points>\n");  
1624    
1625    /* now for the cells */    /* write out the DataArray element for the types */
1626    fprintf(fileHandle_p, "<Cells>\n");    fprintf(fileHandle_p, "<DataArray Name=\"types\" type=\"UInt8\" format=\"ascii\">\n");
1627    /* fprintf(fileHandle_p,    for (i=0; i<numCells; i++) fprintf(fileHandle_p, "%d\n", cellType);
1628        "<DataArray type=\"%s\" Name=\"%s\" format=\"ascii\">\n",);    fprintf(fileHandle_p, "</DataArray>\n");
       */  
   fprintf(fileHandle_p,  
       "</DataArray>\n");  
   fprintf(fileHandle_p, "</Cells>\n");  
1629    
1630    /* finish off the piece */    /* finish off the <Cells> element */
1631    fprintf(fileHandle_p, "</Piece>\n");    fprintf(fileHandle_p, "</Cells>\n");
1632    
1633    /* positions */    /* cell data */
1634    fprintf(fileHandle_p, "object 1 class array type float rank 1 shape %d items %d data follows\n", nDim, mesh_p->Nodes->reducedNumNodes);    if (write_celldata)
1635    for (i = 0; i < mesh_p->Nodes->numNodes; i++) {    {
1636      if (mesh_p->Nodes->toReduced[i]>=0) {      /* mark the active data arrays */
1637         fprintf(fileHandle_p, "%e", mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)]);      bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
1638         for (j = 1; j < nDim; j++) fprintf(fileHandle_p, " %f",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);      fprintf(fileHandle_p, "<CellData");
1639         fprintf(fileHandle_p, "\n");      for (i_data =0 ;i_data<num_data;++i_data)
1640      }      {
1641    }        if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data])
1642    /* connections */        {
1643    /* get a pointer to the relevant mesh component */          /* if the rank == 0:   --> scalar data
1644    if (isEmpty(data_p)) {          * if the rank == 1:   --> vector data
1645        elements=mesh_p->Elements;          * if the rank == 2:   --> tensor data
1646    } else {          */
1647        switch(getFunctionSpaceType(data_p)) {          switch(getDataPointRank(data_pp[i_data]))
1648         case(FINLEY_DEGREES_OF_FREEDOM):          {
1649            nodetype=FINLEY_DEGREES_OF_FREEDOM;          case 0:
1650         case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):            if (! set_scalar)
1651            nodetype=FINLEY_REDUCED_DEGREES_OF_FREEDOM;            {
1652         case(FINLEY_NODES):              fprintf(fileHandle_p," Scalars=\"%s\"",names_p[i_data]);
1653            nodetype=FINLEY_NODES;              set_scalar=TRUE;
1654            isCellCentered=FALSE;            }
1655            elements=mesh_p->Elements;            break;
1656            break;          case 1:
1657         case(FINLEY_ELEMENTS):            if (! set_vector)
1658            isCellCentered=TRUE;            {
1659            elements=mesh_p->Elements;              fprintf(fileHandle_p," Vectors=\"%s\"",names_p[i_data]);
1660            break;              set_vector=TRUE;
1661         case(FINLEY_FACE_ELEMENTS):            }
1662            isCellCentered=TRUE;            break;
1663            elements=mesh_p->FaceElements;          case 2:
1664            break;            if (! set_tensor)
1665         case(FINLEY_POINTS):            {
1666            isCellCentered=TRUE;              fprintf(fileHandle_p," Tensors=\"%s\"",names_p[i_data]);
1667            elements=mesh_p->Points;              set_tensor=TRUE;
1668            break;            }
1669         case(FINLEY_CONTACT_ELEMENTS_1):            break;
1670         case(FINLEY_CONTACT_ELEMENTS_2):          default:
1671            isCellCentered=TRUE;            sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);
1672            elements=mesh_p->ContactElements;            Finley_setError(VALUE_ERROR,error_msg);
1673            break;            fclose(fileHandle_p);
1674         default:            return;
1675            Finley_ErrorCode=TYPE_ERROR;          }
1676            sprintf(Finley_ErrorMsg,"Finley does not know anything about function space type %d",getFunctionSpaceType(data_p));        }
1677            break;      }
1678       }      fprintf(fileHandle_p, ">\n");
1679    }      /* write the arrays */
1680    /* if no element table is present jump over the connection table */      for (i_data =0 ;i_data<num_data;++i_data)
1681    if (elements!=NULL) {      {
1682         ElementTypeId TypeId = elements->ReferenceElement->Type->TypeId;        if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data])
1683         if (TypeId==Line2 || TypeId==Line3 || TypeId==Line4 ) {        {
1684            numVTKNodesPerElement=2;          numPointsPerSample = elements->ReferenceElement->numQuadNodes;
1685            resortIndex=resort[0];          rank = getDataPointRank(data_pp[i_data]);
1686            strcpy(elemTypeStr, "lines");          nComp = getDataPointSize(data_pp[i_data]);
1687         } else if (TypeId==Tri3 || TypeId==Tri6 || TypeId==Tri9 || TypeId==Tri10 ) {          nCompReqd=1;   /* the number of components required by vtk */
1688            numVTKNodesPerElement = 3;          shape=0;
1689            resortIndex=resort[1];          if (rank == 0)
1690            strcpy(elemTypeStr, "triangles");          {
1691         } else if (TypeId==Rec4 || TypeId==Rec8 || TypeId==Rec9 || TypeId==Rec12 || TypeId==Rec16 ) {            nCompReqd = 1;
1692            numVTKNodesPerElement = 4;          }
1693            resortIndex=resort[3];          else if (rank == 1)
1694            strcpy(elemTypeStr, "quads");          {
1695          } else if (TypeId==Tet4 || TypeId==Tet10 || TypeId==Tet16 ) {            shape=getDataPointShape(data_pp[i_data], 0);
1696            numVTKNodesPerElement = 4;            if  (shape>3)
1697            resortIndex=resort[2];            {
1698            strcpy(elemTypeStr, "tetrahedra");              Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");
1699          } else if (TypeId==Hex8 || TypeId==Hex20 || TypeId==Hex32 ) {              fclose(fileHandle_p);
1700            numVTKNodesPerElement = 8;              return;
1701            resortIndex=resort[4];            }
1702            strcpy(elemTypeStr, "cubes");            nCompReqd = 3;
1703          } else {          }
1704            Finley_ErrorCode=VALUE_ERROR;          else
1705            sprintf(Finley_ErrorMsg, "Element type %s is not supported by VTK",elements->ReferenceElement->Type->Name);          {
1706            return;            shape=getDataPointShape(data_pp[i_data], 0);
1707          }            if  (shape>3 || shape != getDataPointShape(data_pp[i_data], 1))
1708          int NN=elements->ReferenceElement->Type->numNodes;            {
1709          fprintf(fileHandle_p, "object 2 class array type int rank 1 shape %d items %d data follows\n",numVTKNodesPerElement, elements->numElements);              Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
1710          for (i = 0; i < elements->numElements; i++) {              fclose(fileHandle_p);
1711            fprintf(fileHandle_p,"%d",mesh_p->Nodes->toReduced[mesh_p->Elements->Nodes[INDEX2(resortIndex[0], i, NN)]]);              return;
1712            for (j = 1; j < numVTKNodesPerElement; j++) {            }
1713               fprintf(fileHandle_p," %d",mesh_p->Nodes->toReduced[mesh_p->Elements->Nodes[INDEX2(resortIndex[j], i, NN)]]);            nCompReqd = 9;
1714            }
1715            fprintf(fileHandle_p, "<DataArray Name=\"%s\" type=\"Float32\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);
1716    
1717            double sampleAvg[nComp];
1718            for (i=0; i<numCells; i++)
1719            {
1720              values = getSampleData(data_pp[i_data], i);
1721              /* averaging over the number of points in the sample */
1722              for (k=0; k<nComp; k++)
1723              {
1724                rtmp = 0.;
1725                for (j=0; j<numPointsPerSample; j++) rtmp += values[INDEX2(k,j,nComp)];
1726                sampleAvg[k] = rtmp/numPointsPerSample;
1727              }
1728              /* if the number of required components is more than the number
1729              * of actual components, pad with zeros
1730              */
1731              /* probably only need to get shape of first element */
1732              /* write the data different ways for scalar, vector and tensor */
1733              if (nCompReqd == 1)
1734              {
1735                fprintf(fileHandle_p, " %e", sampleAvg[0]);
1736              }
1737              else if (nCompReqd == 3)
1738              {
1739                /* write out the data */
1740                for (m=0; m<shape; m++) fprintf(fileHandle_p, " %e", sampleAvg[m]);
1741                for (m=0; m<nCompReqd-shape; m++) fprintf(fileHandle_p, " %e", 0.);
1742              }
1743              else if (nCompReqd == 9)
1744              {
1745                /* tensor data, so have a 3x3 matrix to output as a row
1746                * of 9 data points */
1747                count = 0;
1748                for (m=0; m<shape; m++)
1749                {
1750                  for (n=0; n<shape; n++)
1751                  {
1752                    fprintf(fileHandle_p, " %e", sampleAvg[count]);
1753                    count++;
1754                  }
1755                  for (n=0; n<3-shape; n++) fprintf(fileHandle_p, " %e", 0.);
1756                }
1757                for (m=0; m<3-shape; m++)
1758                  for (n=0; n<3; n++) fprintf(fileHandle_p, " %e", 0.);
1759            }            }
1760            fprintf(fileHandle_p, "\n");            fprintf(fileHandle_p, "\n");
1761          }          }
1762          fprintf(fileHandle_p, "attribute \"element type\" string \"%s\"\n",elemTypeStr);          fprintf(fileHandle_p, "</DataArray>\n");
         fprintf(fileHandle_p, "attribute \"ref\" string \"positions\"\n");  
   
   }  
   /* data */  
   if (!isEmpty(data_p)) {  
       int rank=getDataPointRank(data_p);  
       int* shape=getDataPointShape(data_p);  
       int nComp=getDataPointSize(data_p);  
       fprintf(fileHandle_p, "object 3 class array type float rank %d ", rank);  
       if (0 < rank) {  
          fprintf(fileHandle_p, "shape ");  
          for (i = 0; i < rank; i++) fprintf(fileHandle_p, "%d ", shape[i]);  
       }  
       if (isCellCentered) {  
           int numPointsPerSample=elements->ReferenceElement->numQuadNodes;  
           if (numPointsPerSample) {  
              fprintf(fileHandle_p, "items %d data follows\n", elements->numElements);  
              for (i=0;i<elements->numElements;i++) {  
                  values=getSampleData(data_p,i);  
                  for (k=0;k<nComp;k++) {  
                      rtmp=0.;  
                      for (j=0;j<numPointsPerSample;j++) rtmp+=values[INDEX2(k,j,nComp)];  
                      fprintf(fileHandle_p, " %f", rtmp/numPointsPerSample);  
                  }  
              fprintf(fileHandle_p, "\n");  
              }  
              fprintf(fileHandle_p, "attribute \"dep\" string \"connections\"\n");  
          }  
       } else {  
           fprintf(fileHandle_p, "items %d data follows\n", mesh_p->Nodes->reducedNumNodes);  
           for (i=0;i<mesh_p->Nodes->numNodes;i++) {  
               if (mesh_p->Nodes->toReduced[i]>=0) {  
                  switch (nodetype) {  
                     case(FINLEY_DEGREES_OF_FREEDOM):  
                        values=getSampleData(data_p,mesh_p->Nodes->degreeOfFreedom[i]);  
                        break;  
                     case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):  
                        values=getSampleData(data_p,mesh_p->Nodes->reducedDegreeOfFreedom[i]);  
                        break;  
                     case(FINLEY_NODES):  
                        values=getSampleData(data_p,i);  
                        break;  
                  }  
               }  
               for (k=0;k<nComp;k++) fprintf(fileHandle_p, " %f", values[k]);  
           fprintf(fileHandle_p, "\n");  
           }  
           fprintf(fileHandle_p, "attribute \"dep\" string \"positions\"\n");  
1763        }        }
1764        }
1765        fprintf(fileHandle_p, "</CellData>\n");
1766    }    }
1767      /* point data */
1768    /* and finish it up */    if (write_pointdata)
1769    fprintf(fileHandle_p, "object 4 class field\n");    {
1770    fprintf(fileHandle_p, "component \"positions\" value 1\n");      /* mark the active data arrays */
1771    if (elements!=NULL) fprintf(fileHandle_p, "component \"connections\" value 2\n");      bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
1772    if (!isEmpty(data_p)) fprintf(fileHandle_p, "component \"data\" value 3\n");      fprintf(fileHandle_p, "<PointData");
1773    fprintf(fileHandle_p, "end\n");      for (i_data =0 ;i_data<num_data;++i_data)
1774        {
1775          if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])
1776          {
1777            /* if the rank == 0:   --> scalar data
1778            * if the rank == 1:   --> vector data
1779            * if the rank == 2:   --> tensor data
1780            */
1781            switch(getDataPointRank(data_pp[i_data]))
1782            {
1783            case 0:
1784              if (! set_scalar)
1785              {
1786                fprintf(fileHandle_p," Scalars=\"%s\"",names_p[i_data]);
1787                set_scalar=TRUE;
1788              }
1789              break;
1790            case 1:
1791              if (! set_vector)
1792              {
1793                fprintf(fileHandle_p," Vectors=\"%s\"",names_p[i_data]);
1794                set_vector=TRUE;
1795              }
1796              break;
1797            case 2:
1798              if (! set_tensor)
1799              {
1800                fprintf(fileHandle_p," Tensors=\"%s\"",names_p[i_data]);
1801                set_tensor=TRUE;
1802              }
1803              break;
1804            default:
1805              sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);
1806              Finley_setError(VALUE_ERROR,error_msg);
1807              fclose(fileHandle_p);
1808              return;
1809            }
1810          }
1811        }
1812        fprintf(fileHandle_p, ">\n");
1813        /* write the arrays */
1814        for (i_data =0 ;i_data<num_data;++i_data)
1815        {
1816          if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])
1817          {
1818            numPointsPerSample = elements->ReferenceElement->numQuadNodes;
1819            rank = getDataPointRank(data_pp[i_data]);
1820            nComp = getDataPointSize(data_pp[i_data]);
1821            nCompReqd=1;   /* the number of components required by vtk */
1822            shape=0;
1823            if (rank == 0)
1824            {
1825              nCompReqd = 1;
1826            }
1827            else if (rank == 1)
1828            {
1829              shape=getDataPointShape(data_pp[i_data], 0);
1830              if  (shape>3)
1831              {
1832                Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");
1833                fclose(fileHandle_p);
1834                return;
1835              }
1836              nCompReqd = 3;
1837            }
1838            else
1839            {
1840              shape=getDataPointShape(data_pp[i_data], 0);
1841              if  (shape>3 || shape != getDataPointShape(data_pp[i_data], 1))
1842              {
1843                Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
1844                fclose(fileHandle_p);
1845                return;
1846              }
1847              nCompReqd = 9;
1848            }
1849            fprintf(fileHandle_p, "<DataArray Name=\"%s\" type=\"Float32\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);
1850            /* write out the data */
1851            /* if the number of required components is more than the number
1852            * of actual components, pad with zeros
1853            */
1854            bool_t do_write=TRUE;
1855            for (i=0; i<mesh_p->Nodes->numNodes; i++)
1856            {
1857              if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1858              {
1859                if (mesh_p->Nodes->toReduced[i]>=0)
1860                {
1861                  switch(getFunctionSpaceType(data_pp[i_data]))
1862                  {
1863                  case FINLEY_DEGREES_OF_FREEDOM:
1864                    values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);
1865                    break;
1866                  case FINLEY_REDUCED_DEGREES_OF_FREEDOM:
1867                    values = getSampleData(data_pp[i_data],mesh_p->Nodes->reducedDegreeOfFreedom[i]);
1868                    break;
1869                  case FINLEY_NODES:
1870                    values = getSampleData(data_pp[i_data],i);
1871                    break;
1872                  }
1873                  do_write=TRUE;
1874                }
1875                else
1876                {
1877                  do_write=FALSE;
1878                }
1879              }
1880              else
1881              {
1882                do_write=TRUE;
1883                switch(getFunctionSpaceType(data_pp[i_data]))
1884                {
1885                case FINLEY_DEGREES_OF_FREEDOM:
1886                  values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);
1887                  break;
1888                case FINLEY_NODES:
1889                  values = getSampleData(data_pp[i_data],i);
1890                  break;
1891                }
1892              }
1893              /* write the data different ways for scalar, vector and tensor */
1894              if (do_write)
1895              {
1896                if (nCompReqd == 1)
1897                {
1898                  fprintf(fileHandle_p, " %e", values[0]);
1899                }
1900                else if (nCompReqd == 3)
1901                {
1902                  for (m=0; m<shape; m++) fprintf(fileHandle_p, " %e", values[m]);
1903                  for (m=0; m<nCompReqd-shape; m++) fprintf(fileHandle_p, " %e", 0.);
1904                }
1905                else if (nCompReqd == 9)
1906                {
1907                  /* tensor data, so have a 3x3 matrix to output as a row
1908                  * of 9 data points */
1909                  count = 0;
1910                  for (m=0; m<shape; m++)
1911                  {
1912                    for (n=0; n<shape; n++)
1913                    {
1914                      fprintf(fileHandle_p, " %e", values[count]);
1915                      count++;
1916                    }
1917                    for (n=0; n<3-shape; n++) fprintf(fileHandle_p, " %e", 0.);
1918                  }
1919                  for (m=0; m<3-shape; m++)
1920                    for (n=0; n<3; n++) fprintf(fileHandle_p, " %e", 0.);
1921                }
1922                fprintf(fileHandle_p, "\n");
1923              }
1924            }
1925            fprintf(fileHandle_p, "</DataArray>\n");
1926          }
1927        }
1928        fprintf(fileHandle_p, "</PointData>\n");
1929      }
1930      /* finish off the piece */
1931      fprintf(fileHandle_p, "</Piece>\n");
1932    
1933    fprintf(fileHandle_p, "</UnstructuredGrid>\n");    fprintf(fileHandle_p, "</UnstructuredGrid>\n");
1934    /* write the xml footer */    /* write the xml footer */
# Line 234  void Finley_Mesh_saveVTK(const char * fi Line 1937  void Finley_Mesh_saveVTK(const char * fi
1937    fclose(fileHandle_p);    fclose(fileHandle_p);
1938    return;    return;
1939  }  }
1940    #endif
1941    
 /*  
  * $Log$  
  * Revision 1.2  2005/02/14 04:14:42  jgs  
  * *** empty log message ***  
  *  
  * 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.110  
changed lines
  Added in v.818

  ViewVC Help
Powered by ViewVC 1.1.26