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

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

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

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

Legend:
Removed from v.155  
changed lines
  Added in v.794

  ViewVC Help
Powered by ViewVC 1.1.26