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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26