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

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

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

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