/[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 794 - (hide annotations)
Sun Jul 30 03:45:01 2006 UTC (13 years, 2 months ago) by dhawcroft
File MIME type: text/plain
File size: 60177 byte(s)


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 794 /* MPI-IO version: Derick Hawcroft, d.hawcroft@uq.edu.au */
22 dhawcroft 793
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 dhawcroft 794 MPI version notes:
33 dhawcroft 793
34 dhawcroft 794 ******************************************************************************
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 dhawcroft 793 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 dhawcroft 794
46 dhawcroft 793 */
47 jgs 150
48 dhawcroft 793 #ifdef PASO_MPI
49    
50 dhawcroft 794
51 dhawcroft 793 //#define MPIO_HINTS
52    
53    
54 dhawcroft 794
55 dhawcroft 793 #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 dhawcroft 794 numInternalCells,
84     numBoundaryCells;
85 dhawcroft 793
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 jgs 150 char error_msg[LenErrorMsg_MAX];
100 dhawcroft 793
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 dhawcroft 794 numLocalNodes,
119     numBoundaryNodes,
120     localDOF;
121    
122 dhawcroft 793
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 dhawcroft 794 struct localIndexCache
143 dhawcroft 793 {
144     index_t *values;
145     int size;
146     };
147 dhawcroft 794 typedef struct localIndexCache localIndexCache;
148 dhawcroft 793
149 dhawcroft 794 localIndexCache nodeCache,
150     elementCache;
151 dhawcroft 793
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 dhawcroft 794 numBoundaryNodes = dist->numBoundary;
302 dhawcroft 793
303 dhawcroft 794 localDOF = dist->numLocal;
304 dhawcroft 793
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 dhawcroft 794 numLocalNodes=localDOF;
466    
467 dhawcroft 793 // 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 dhawcroft 794 DOFNodes = MEMALLOC(mesh_p->Nodes->numNodes,index_t);
478     nodeCache.values = MEMALLOC( numLocalNodes, index_t);
479     index_t bc_pos = 0;
480 dhawcroft 793 for (i = 0; i < mesh_p->Nodes->numNodes; i++)
481 dhawcroft 794
482 dhawcroft 793 {
483 dhawcroft 794 // 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 dhawcroft 793 DOFNodes[mesh_p->Nodes->degreeOfFreedom[i]] = i;
486 dhawcroft 794
487 dhawcroft 793 /* 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 dhawcroft 794
512     nodesGlobal = MEMALLOC(mesh_p->Nodes->numNodes ,index_t);
513 dhawcroft 793
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 dhawcroft 794
581     // TODO: Improve on upper bound
582     size_t sz = numLocalCells*6*numVTKNodesPerElement + numLocalCells;
583     char *cellBuf = MEMALLOC(sz,char);
584 dhawcroft 793 cellBuf[0] = '\0';
585     tsz=0;
586     pos = 0;
587 dhawcroft 794 // numCells?
588     elementCache.values = MEMALLOC(numLocalCells,index_t);
589 dhawcroft 793 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 dhawcroft 794 char tmpbuf2[20*20];
611 dhawcroft 793 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 dhawcroft 794
678 dhawcroft 793 MPI_File_write_ordered(fh, cellBuf,tsz, MPI_CHAR, &status);
679     MEMFREE(cellBuf);
680     MPIO_DEBUG(" Done Writing Connectivity ")
681 dhawcroft 794 MPIO_DEBUG(" Writing Offsets & Types... ")
682 dhawcroft 793
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 dhawcroft 794 sz += numGlobalCells;
698 dhawcroft 793
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 dhawcroft 794 MPIO_DEBUG(" Done Writing Offsets & Types ")
731 dhawcroft 793
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 dhawcroft 794 // closing celldata tag
1158 dhawcroft 793 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];
1195 jgs 110 /* if there is no mesh we just return */
1196     if (mesh_p==NULL) return;
1197 jgs 153
1198 gross 687 int i, j, k, numVTKNodesPerElement,i_data,m, count, n, rank,shape, numPoints, cellType, numCells,
1199 dhawcroft 793 nDim, numPointsPerSample, nComp, nCompReqd;
1200 gross 687
1201 jgs 150 index_t j2;
1202 jgs 113 double* values, rtmp;
1203 gross 687 char elemTypeStr[32];
1204 dhawcroft 793
1205 jgs 153 /* open the file and check handle */
1206    
1207     FILE * fileHandle_p = fopen(filename_p, "w");
1208 dhawcroft 793 if (fileHandle_p==NULL)
1209     {
1210 jgs 153 sprintf(error_msg, "saveVTK: File %s could not be opened for writing.", filename_p);
1211     Finley_setError(IO_ERROR,error_msg);
1212     return;
1213 jgs 150 }
1214 jgs 153 /* find the mesh type to be written */
1215     int nodetype=FINLEY_DEGREES_OF_FREEDOM;
1216     int elementtype=FINLEY_UNKNOWN;
1217     bool_t isCellCentered[num_data],write_celldata=FALSE,write_pointdata=FALSE;
1218 dhawcroft 793 for (i_data=0;i_data<num_data;++i_data)
1219     {
1220     if (! isEmpty(data_pp[i_data]))
1221     {
1222     switch(getFunctionSpaceType(data_pp[i_data]))
1223     {
1224     case FINLEY_DEGREES_OF_FREEDOM:
1225     nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1226     if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)
1227     {
1228     elementtype=FINLEY_ELEMENTS;
1229 jgs 153 }
1230 dhawcroft 793 else
1231     {
1232     Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1233     fclose(fileHandle_p);
1234     return;
1235 jgs 153 }
1236 dhawcroft 793 isCellCentered[i_data]=FALSE;
1237     break;
1238     case FINLEY_REDUCED_DEGREES_OF_FREEDOM:
1239     nodetype = FINLEY_REDUCED_DEGREES_OF_FREEDOM;
1240     if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)
1241     {
1242     elementtype=FINLEY_ELEMENTS;
1243     }
1244     else
1245     {
1246     Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1247     fclose(fileHandle_p);
1248     return;
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 jgs 153 }
1352     /* select nomber of points and the mesh component */
1353 gross 687 numPoints = mesh_p->Nodes->numNodes;
1354 dhawcroft 793 if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1355     {
1356     numPoints = mesh_p->Nodes->reducedNumNodes;
1357 jgs 153 }
1358 dhawcroft 793 else
1359     {
1360     numPoints = mesh_p->Nodes->numNodes;
1361     }
1362 jgs 153 if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;
1363     Finley_ElementFile* elements=NULL;
1364 dhawcroft 793 switch(elementtype)
1365     {
1366     case FINLEY_ELEMENTS:
1367     elements=mesh_p->Elements;
1368     break;
1369     case FINLEY_FACE_ELEMENTS:
1370     elements=mesh_p->FaceElements;
1371     break;
1372     case FINLEY_POINTS:
1373     elements=mesh_p->Points;
1374     break;
1375     case FINLEY_CONTACT_ELEMENTS_1:
1376     elements=mesh_p->ContactElements;
1377     break;
1378 jgs 113 }
1379 dhawcroft 793 if (elements==NULL)
1380     {
1381     Finley_setError(SYSTEM_ERROR,"saveVTK: undefined element file");
1382     fclose(fileHandle_p);
1383     return;
1384 jgs 113 }
1385 jgs 153 /* map finley element type to VTK element type */
1386 dhawcroft 793 numCells = elements->numElements;
1387 jgs 153 ElementTypeId TypeId;
1388 dhawcroft 793 if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1389     {
1390     TypeId = elements->LinearReferenceElement->Type->TypeId;
1391 jgs 113 }
1392 dhawcroft 793 else
1393     {
1394     TypeId = elements->ReferenceElement->Type->TypeId;
1395     }
1396 jgs 110
1397 dhawcroft 793 switch(TypeId)
1398     {
1399     case Point1:
1400     case Line2Face:
1401     case Line3Face:
1402     case Point1_Contact:
1403     case Line2Face_Contact:
1404     case Line3Face_Contact:
1405     cellType = VTK_VERTEX;
1406     numVTKNodesPerElement = 1;
1407     strcpy(elemTypeStr, "VTK_VERTEX");
1408     break;
1409 jgs 153
1410 dhawcroft 793 case Line2:
1411     case Tri3Face:
1412     case Rec4Face:
1413     case Line2_Contact:
1414     case Tri3_Contact:
1415     case Tri3Face_Contact:
1416     case Rec4Face_Contact:
1417     cellType = VTK_LINE;
1418     numVTKNodesPerElement = 2;
1419     strcpy(elemTypeStr, "VTK_LINE");
1420     break;
1421 jgs 153
1422 dhawcroft 793 case Tri3:
1423     case Tet4Face:
1424     case Tet4Face_Contact:
1425     cellType = VTK_TRIANGLE;
1426     numVTKNodesPerElement = 3;
1427     strcpy(elemTypeStr, "VTK_TRIANGLE");
1428     break;
1429 jgs 153
1430 dhawcroft 793 case Rec4:
1431     case Hex8Face:
1432     case Rec4_Contact:
1433     case Hex8Face_Contact:
1434     cellType = VTK_QUAD;
1435     numVTKNodesPerElement = 4;
1436     strcpy(elemTypeStr, "VTK_QUAD");
1437     break;
1438 jgs 153
1439 dhawcroft 793 case Tet4:
1440     cellType = VTK_TETRA;
1441     numVTKNodesPerElement = 4;
1442     strcpy(elemTypeStr, "VTK_TETRA");
1443     break;
1444 jgs 153
1445 dhawcroft 793 case Hex8:
1446     cellType = VTK_HEXAHEDRON;
1447     numVTKNodesPerElement = 8;
1448     strcpy(elemTypeStr, "VTK_HEXAHEDRON");
1449     break;
1450 jgs 153
1451 dhawcroft 793 case Line3:
1452     case Tri6Face:
1453     case Rec8Face:
1454     case Line3_Contact:
1455     case Tri6Face_Contact:
1456     case Rec8Face_Contact:
1457     cellType = VTK_QUADRATIC_EDGE;
1458     numVTKNodesPerElement = 3;
1459     strcpy(elemTypeStr, "VTK_QUADRATIC_EDGE");
1460     break;
1461 jgs 153
1462 dhawcroft 793 case Tri6:
1463     case Tet10Face:
1464     case Tri6_Contact:
1465     case Tet10Face_Contact:
1466     cellType = VTK_QUADRATIC_TRIANGLE;
1467     numVTKNodesPerElement = 6;
1468     strcpy(elemTypeStr, "VTK_QUADRATIC_TRIANGLE");
1469     break;
1470 jgs 153
1471 dhawcroft 793 case Rec8:
1472     case Hex20Face:
1473     case Rec8_Contact:
1474     case Hex20Face_Contact:
1475     cellType = VTK_QUADRATIC_QUAD;
1476     numVTKNodesPerElement = 8;
1477     strcpy(elemTypeStr, "VTK_QUADRATIC_QUAD");
1478     break;
1479 jgs 153
1480 dhawcroft 793 case Tet10:
1481     cellType = VTK_QUADRATIC_TETRA;
1482     numVTKNodesPerElement = 10;
1483     strcpy(elemTypeStr, "VTK_QUADRATIC_TETRA");
1484     break;
1485 jgs 153
1486 dhawcroft 793 case Hex20:
1487     cellType = VTK_QUADRATIC_HEXAHEDRON;
1488     numVTKNodesPerElement = 20;
1489     strcpy(elemTypeStr, "VTK_QUADRATIC_HEXAHEDRON");
1490     break;
1491 jgs 153
1492 dhawcroft 793 default:
1493     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 jgs 153 /* xml header */
1499     fprintf(fileHandle_p, "<?xml version=\"1.0\"?>\n");
1500     fprintf(fileHandle_p, "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\">\n");
1501 jgs 110
1502 jgs 153 /* finley uses an unstructured mesh, so UnstructuredGrid *should* work */
1503     fprintf(fileHandle_p, "<UnstructuredGrid>\n");
1504    
1505     /* is there only one "piece" to the data?? */
1506     fprintf(fileHandle_p, "<Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n",numPoints, numCells);
1507     /* now for the points; equivalent to positions section in saveDX() */
1508     /* "The points element explicitly defines coordinates for each point
1509 dhawcroft 793 * individually. It contains one DataArray element describing an array
1510     * with three components per value, each specifying the coordinates of one
1511     * point" - from Vtk User's Guide
1512     */
1513 jgs 153 fprintf(fileHandle_p, "<Points>\n");
1514 dhawcroft 793 /*
1515     * the reason for this if statement is explained in the long comment below
1516     */
1517 gross 687 nDim = mesh_p->Nodes->numDim;
1518 jgs 153 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
1520 dhawcroft 793 * 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
1522     * dimension, and keep the visualisers happy.
1523     * Indeed, if nDim is less than 3, must pad all empty dimensions, so
1524     * that the total number of dims is 3.
1525 jgs 153 */
1526 dhawcroft 793 if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1527     {
1528     for (i = 0; i < mesh_p->Nodes->numNodes; i++)
1529     {
1530     if (mesh_p->Nodes->toReduced[i]>=0)
1531     {
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     fprintf(fileHandle_p, "\n");
1535     }
1536     }
1537     }
1538     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 jgs 153 fprintf(fileHandle_p, "</DataArray>\n");
1550     fprintf(fileHandle_p, "</Points>\n");
1551 jgs 113
1552 jgs 153 /* write out the DataArray element for the connectivity */
1553 jgs 113
1554 jgs 153 int NN = elements->ReferenceElement->Type->numNodes;
1555     fprintf(fileHandle_p, "<Cells>\n");
1556     fprintf(fileHandle_p, "<DataArray Name=\"connectivity\" type=\"Int32\" format=\"ascii\">\n");
1557 jgs 113
1558 dhawcroft 793 if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1559     {
1560     for (i = 0; i < numCells; i++)
1561     {
1562     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     fprintf(fileHandle_p, "\n");
1565     }
1566     }
1567     else if (VTK_QUADRATIC_HEXAHEDRON==cellType)
1568     {
1569     for (i = 0; i < numCells; i++)
1570     {
1571     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(0, i, NN)],
1573     elements->Nodes[INDEX2(1, i, NN)],
1574     elements->Nodes[INDEX2(2, i, NN)],
1575     elements->Nodes[INDEX2(3, i, NN)],
1576     elements->Nodes[INDEX2(4, i, NN)],
1577     elements->Nodes[INDEX2(5, i, NN)],
1578     elements->Nodes[INDEX2(6, i, NN)],
1579     elements->Nodes[INDEX2(7, i, NN)],
1580     elements->Nodes[INDEX2(8, i, NN)],
1581     elements->Nodes[INDEX2(9, i, NN)],
1582     elements->Nodes[INDEX2(10, i, NN)],
1583     elements->Nodes[INDEX2(11, i, NN)],
1584     elements->Nodes[INDEX2(16, i, NN)],
1585     elements->Nodes[INDEX2(17, i, NN)],
1586     elements->Nodes[INDEX2(18, i, NN)],
1587     elements->Nodes[INDEX2(19, i, NN)],
1588     elements->Nodes[INDEX2(12, i, NN)],
1589     elements->Nodes[INDEX2(13, i, NN)],
1590     elements->Nodes[INDEX2(14, i, NN)],
1591     elements->Nodes[INDEX2(15, i, NN)]);
1592     }
1593     }
1594     else if (numVTKNodesPerElement!=NN)
1595     {
1596     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 jgs 153 fprintf(fileHandle_p, "</DataArray>\n");
1611 jgs 110
1612 jgs 153 /* write out the DataArray element for the offsets */
1613     fprintf(fileHandle_p, "<DataArray Name=\"offsets\" type=\"Int32\" format=\"ascii\">\n");
1614     for (i=numVTKNodesPerElement; i<=numCells*numVTKNodesPerElement; i+=numVTKNodesPerElement) fprintf(fileHandle_p, "%d\n", i);
1615     fprintf(fileHandle_p, "</DataArray>\n");
1616    
1617     /* write out the DataArray element for the types */
1618     fprintf(fileHandle_p, "<DataArray Name=\"types\" type=\"UInt8\" format=\"ascii\">\n");
1619     for (i=0; i<numCells; i++) fprintf(fileHandle_p, "%d\n", cellType);
1620     fprintf(fileHandle_p, "</DataArray>\n");
1621    
1622     /* finish off the <Cells> element */
1623     fprintf(fileHandle_p, "</Cells>\n");
1624    
1625     /* cell data */
1626 dhawcroft 793 if (write_celldata)
1627     {
1628     /* mark the active data arrays */
1629     bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
1630     fprintf(fileHandle_p, "<CellData");
1631     for (i_data =0 ;i_data<num_data;++i_data)
1632     {
1633     if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data])
1634     {
1635     /* if the rank == 0: --> scalar data
1636     * if the rank == 1: --> vector data
1637     * if the rank == 2: --> tensor data
1638     */
1639     switch(getDataPointRank(data_pp[i_data]))
1640     {
1641     case 0:
1642     if (! set_scalar)
1643     {
1644     fprintf(fileHandle_p," Scalars=\"%s\"",names_p[i_data]);
1645     set_scalar=TRUE;
1646     }
1647     break;
1648     case 1:
1649     if (! set_vector)
1650     {
1651     fprintf(fileHandle_p," Vectors=\"%s\"",names_p[i_data]);
1652     set_vector=TRUE;
1653     }
1654     break;
1655     case 2:
1656     if (! set_tensor)
1657     {
1658     fprintf(fileHandle_p," Tensors=\"%s\"",names_p[i_data]);
1659     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 jgs 150 }
1749 dhawcroft 793 for (m=0; m<3-shape; m++)
1750     for (n=0; n<3; n++) fprintf(fileHandle_p, " %e", 0.);
1751     }
1752     fprintf(fileHandle_p, "\n");
1753     }
1754     fprintf(fileHandle_p, "</DataArray>\n");
1755     }
1756     }
1757     fprintf(fileHandle_p, "</CellData>\n");
1758 jgs 110 }
1759 jgs 153 /* point data */
1760 dhawcroft 793 if (write_pointdata)
1761     {
1762     /* mark the active data arrays */
1763     bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
1764     fprintf(fileHandle_p, "<PointData");
1765     for (i_data =0 ;i_data<num_data;++i_data)
1766     {
1767     if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])
1768     {
1769     /* if the rank == 0: --> scalar data
1770     * if the rank == 1: --> vector data
1771     * if the rank == 2: --> tensor data
1772     */
1773     switch(getDataPointRank(data_pp[i_data]))
1774     {
1775     case 0:
1776     if (! set_scalar)
1777     {
1778     fprintf(fileHandle_p," Scalars=\"%s\"",names_p[i_data]);
1779     set_scalar=TRUE;
1780     }
1781     break;
1782     case 1:
1783     if (! set_vector)
1784     {
1785     fprintf(fileHandle_p," Vectors=\"%s\"",names_p[i_data]);
1786     set_vector=TRUE;
1787     }
1788     break;
1789     case 2:
1790     if (! set_tensor)
1791     {
1792     fprintf(fileHandle_p," Tensors=\"%s\"",names_p[i_data]);
1793     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 jgs 153 }
1909 dhawcroft 793 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 jgs 153 }
1914 dhawcroft 793 fprintf(fileHandle_p, "\n");
1915     }
1916     }
1917     fprintf(fileHandle_p, "</DataArray>\n");
1918     }
1919     }
1920     fprintf(fileHandle_p, "</PointData>\n");
1921 jgs 153 }
1922 jgs 113 /* finish off the piece */
1923     fprintf(fileHandle_p, "</Piece>\n");
1924 jgs 110
1925     fprintf(fileHandle_p, "</UnstructuredGrid>\n");
1926     /* write the xml footer */
1927     fprintf(fileHandle_p, "</VTKFile>\n");
1928     /* close the file */
1929     fclose(fileHandle_p);
1930     return;
1931     }
1932 dhawcroft 793 #endif
1933    

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26