/[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 1030 - (hide annotations)
Wed Mar 14 05:14:44 2007 UTC (12 years, 7 months ago) by phornby
File MIME type: text/plain
File size: 61733 byte(s)
Assemble_CopyElementData.c - resolve with my local edits.
Mesh_saveVTK.c - remove debug print.
DataArrayView.h - remove dllim/export declarations on template functions.


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 796
46 dhawcroft 818 Notable Notables:
47     the struct localIndexCache: stores local domain indices for faster reference
48 dhawcroft 793 */
49 jgs 150
50 dhawcroft 793 #ifdef PASO_MPI
51    
52 dhawcroft 794
53 dhawcroft 793 //#define MPIO_HINTS
54    
55    
56 dhawcroft 794
57 dhawcroft 793 #define MPIO_DEBUG(str) \
58     { \
59     if(myRank == 0) \
60     printf("==== MPI-IO => %s \n", str); \
61     }
62    
63     void Finley_Mesh_saveVTK_MPIO(const char * filename_p, Finley_Mesh *mesh_p, const dim_t num_data,char* *names_p, escriptDataC* *data_pp)
64     {
65     int numPoints,
66     numCells = -1,
67     myRank,comm,gsize,
68     numLocal,
69     nDim,
70     shape;
71 dhawcroft 796 size_t __n;
72 dhawcroft 793 int i,j,k,m,n,count;
73     int numGlobalCells = 0;
74     index_t *nodesGlobal=NULL; // used to get the connectivity right for VTK
75    
76     /* variables associatted with write_celldata/pointdata */
77     int numPointsPerSample,
78     nComp,
79     nCompReqd;
80     double* values, rtmp;
81    
82     // Local element info (for debugging)
83     size_t numLocalCells,
84 dhawcroft 796 numInternalCells,
85     numBoundaryCells;
86 dhawcroft 793
87     int rank;
88    
89     int amode = MPI_MODE_CREATE | MPI_MODE_WRONLY | MPI_MODE_SEQUENTIAL;
90    
91     comm = mesh_p->Nodes->MPIInfo->comm;
92     myRank = mesh_p->Nodes->MPIInfo->rank;
93     gsize = mesh_p->Nodes->MPIInfo->size;
94    
95     MPI_File fh;
96     MPI_Status status;
97     MPI_Request req;
98     MPI_Info infoHints;
99    
100 jgs 150 char error_msg[LenErrorMsg_MAX];
101 dhawcroft 793
102     int i_data;
103    
104     int nodetype=FINLEY_DEGREES_OF_FREEDOM;
105     int elementtype=FINLEY_UNKNOWN;
106     bool_t isCellCentered[num_data],write_celldata=FALSE,write_pointdata=FALSE;
107    
108     ElementTypeId TypeId;
109    
110     int numVTKNodesPerElement;
111     int cellType;
112     char elemTypeStr[32];
113    
114     Finley_ElementFile* elements=NULL;
115    
116    
117     // Local node info
118     int numInternalNodes,
119 dhawcroft 796 numLocalNodes,
120     numBoundaryNodes,
121 dhawcroft 818 localDOF; // equals to (DOF of Internal Nodes) + (DOF of Boundary Nodes) of local domain
122 dhawcroft 793
123 dhawcroft 796
124 dhawcroft 793 nDim = mesh_p->Nodes->numDim;
125    
126     #ifdef MPIO_HINTS
127     MPI_Info_create(&infoHints);
128     // MPI_Info_set(infoHints, "striping_unit", "424288");
129     // MPI_Info_set(infoHints, "striping_factor", "16");
130     // MPI_Info_set(infoHints, "collective_buffering", "true");
131     // MPI_Info_set(infoHints, "cb_block_size", "131072");
132     // MPI_Info_set(infoHints, "cb_buffer_size", "1048567");
133     // MPI_Info_set(infoHints, "cb_nodes", "8");
134     // MPI_Info_set(infoHints, "access_style", "write_once, sequential");
135    
136     //XFS only
137     // MPI_Info_set(infoHints, "direct_write", "true");
138     #else
139     infoHints = MPI_INFO_NULL;
140     #endif
141    
142 dhawcroft 818 // Holds a local node/element index into the global array
143 dhawcroft 794 struct localIndexCache
144 dhawcroft 793 {
145     index_t *values;
146     int size;
147     };
148    
149 dhawcroft 818 struct 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 gross 850 "<DataArray NumberOfComponents=\"%d\" type=\"Float64\" format=\"ascii\">\n"
456 dhawcroft 793 ,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 dhawcroft 796
467 dhawcroft 818 // we will be writing values which vary from 13-15 chars hence the strlen()
468     char* largebuf = MEMALLOC( numLocalNodes*15*nDim + numLocalNodes*2 + 1 ,char);
469 dhawcroft 793 largebuf[0] = '\0';
470 dhawcroft 818 char tmpbuf[15];
471 dhawcroft 793 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 dhawcroft 796 nodeCache.values = MEMALLOC( numLocalNodes, index_t);
479 dhawcroft 794 index_t bc_pos = 0;
480 dhawcroft 796
481 dhawcroft 818 // Custom string concat: avoids expensive strlen(3) call by strcat(3)
482     // Note the implicit assumption on the variable "tsz"
483 dhawcroft 796 int __len,__j;
484     char *zero = "0.000000e+00 ";
485     char *newline = "\n";
486 dhawcroft 818
487     #define __STRCAT(dest,chunk,tsz) \
488 dhawcroft 796 { \
489 dhawcroft 818 __len = strlen(chunk); \
490 dhawcroft 796 __j = -1; \
491     while(__j++ < __len) \
492 dhawcroft 818 *(dest+tsz+__j)=*(chunk+__j); \
493     tsz+=__len; \
494 dhawcroft 796 }
495 dhawcroft 818
496     // Loop over all nodes
497 dhawcroft 793 for (i = 0; i < mesh_p->Nodes->numNodes; i++)
498     {
499 dhawcroft 794 // This is the bit that will break for periodic BCs because it assumes that there is a one to one
500     // correspondance between nodes and Degrees of freedom
501 dhawcroft 818 //TODO: handle periodic BC's
502 dhawcroft 793 DOFNodes[mesh_p->Nodes->degreeOfFreedom[i]] = i;
503 dhawcroft 796
504 dhawcroft 818 // Is this node local to the domain ?
505 dhawcroft 793 if( mesh_p->Nodes->degreeOfFreedom[i] < localDOF )
506     {
507     for (j = 0; j < nDim; j++)
508     {
509     sprintf(tmpbuf,"%e ", mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)] );
510 dhawcroft 818 __STRCAT(largebuf,tmpbuf,tsz)
511 dhawcroft 793 }
512     for (k=0; k<3-nDim; k++)
513     {
514 dhawcroft 818 __STRCAT(largebuf,zero,tsz)
515 dhawcroft 793 }
516 dhawcroft 818 __STRCAT(largebuf,newline,tsz)
517 dhawcroft 793 nodeCache.values[numNodesOutput++]=i;
518     }
519     }
520    
521     nodeCache.size=numNodesOutput;
522    
523 dhawcroft 796 largebuf[tsz] = '\0';
524 dhawcroft 793 MPI_File_write_ordered(fh, largebuf,tsz, MPI_CHAR, &status);
525     MEMFREE(largebuf);
526 dhawcroft 796
527 dhawcroft 794 nodesGlobal = MEMALLOC(mesh_p->Nodes->numNodes ,index_t);
528 dhawcroft 793
529     // form distribution info on who output which nodes
530     vtxdist = MEMALLOC( gsize+1, index_t );
531     vtxdist[0]=0;
532     MPI_Allgather(&numNodesOutput,1,MPI_INT,vtxdist+1,1,MPI_INT,comm);
533     for( i=0; i<gsize; i++ )
534     vtxdist[i+1]+=vtxdist[i];
535    
536     // will not work for periodic boundary conditions
537     // calculate the local nodes file positions
538     pos = 0;
539     for( i=0; i<mesh_p->Nodes->numNodes; i++ )
540     {
541     if( mesh_p->Nodes->degreeOfFreedom[i]< localDOF )
542     {
543     nodesGlobal[i] = vtxdist[myRank] + pos++;
544     }
545     else
546     nodesGlobal[i] = -1;
547     }
548    
549     // communicate the local Nodes file position to the interested parties
550     // send local info
551     forwardBuffer = MEMALLOC( mesh_p->Nodes->numNodes, index_t );
552     for( n=0; n < dist->numNeighbours; n++ )
553     {
554     if( dist->edges[n]->numForward)
555     {
556     for( i=0; i < dist->edges[n]->numForward; i++ )
557     forwardBuffer[i] = nodesGlobal[DOFNodes[dist->edges[n]->indexForward[i] ]];
558     Paso_CommBuffer_pack( mesh_p->Nodes->CommBuffer, dist->neighbours[n], NULL, forwardBuffer, sizeof(index_t), 0 );
559     Paso_CommBuffer_send( mesh_p->Nodes->CommBuffer, dist->neighbours[n], sizeof(index_t) );
560     }
561     }
562     // receive external info
563     backwardBuffer = MEMALLOC( mesh_p->Nodes->numNodes, index_t );
564     for( n=0; n < dist->numNeighbours; n++ )
565     {
566     if( dist->edges[n]->numBackward )
567     {
568     Paso_CommBuffer_recv(mesh_p->Nodes->CommBuffer, dist->neighbours[n], sizeof(index_t));
569     Paso_CommBuffer_unpack(mesh_p->Nodes->CommBuffer, dist->neighbours[n], NULL, backwardBuffer, sizeof(index_t), 0 );
570 dhawcroft 818 /* TODO: voodoo to handle perdiodic BC's */
571 dhawcroft 793 for( i=0; i<dist->edges[n]->numBackward; i++ )
572     nodesGlobal[DOFNodes[dist->edges[n]->indexBackward[i] ]] = backwardBuffer[i];
573     }
574     }
575 dhawcroft 818
576 dhawcroft 793
577 dhawcroft 818
578 dhawcroft 793 MEMFREE(vtxdist);
579     MEMFREE(DOFNodes);
580     MEMFREE(backwardBuffer);
581     MEMFREE(forwardBuffer);
582    
583     if( myRank == 0)
584     {
585     char* tags = "</DataArray>\n</Points>\n<Cells>\n<DataArray Name=\"connectivity\" type=\"Int32\" " \
586     "format=\"ascii\">\n" ;
587     MPI_File_iwrite_shared(fh,tags,strlen(tags),MPI_CHAR,&req);
588     MPI_Wait(&req,&status);
589     }
590     MPIO_DEBUG(" Done Writing Coordinate Points ")
591    
592     /* BEGIN CONNECTIVITY */
593    
594     int NN = elements->ReferenceElement->Type->numNodes; /* num Nodes holding ref-element */
595    
596     // Collective
597     MPIO_DEBUG(" Writing Connectivity... ")
598 dhawcroft 796
599 dhawcroft 794 size_t sz = numLocalCells*6*numVTKNodesPerElement + numLocalCells;
600 dhawcroft 796 largebuf = MEMALLOC(sz,char);
601     largebuf[0] = '\0';
602 dhawcroft 793 tsz=0;
603     pos = 0;
604 dhawcroft 794 // numCells?
605     elementCache.values = MEMALLOC(numLocalCells,index_t);
606 dhawcroft 793 if (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM)
607     {
608     for (i = 0; i < numCells; i++)
609     {
610 dhawcroft 796
611 dhawcroft 793 if (elements->Id[i] >= elements->elementDistribution->vtxdist[i] && elements->Id[i] <= elements->elementDistribution->vtxdist[i+1] - 1 )
612     {
613     for (j = 0; j < numVTKNodesPerElement; j++)
614     {
615     sprintf(tmpbuf,"%d ",nodesGlobal[mesh_p->Nodes->toReduced[elements->Nodes[INDEX2(elements->ReferenceElement->Type->linearNodes[j], i, NN)]]]);
616 dhawcroft 818 __STRCAT(largebuf,tmpbuf,tsz)
617 dhawcroft 793 }
618 dhawcroft 818 __STRCAT(largebuf,newline,tsz)
619 dhawcroft 793 elementCache.values[pos++]=i;
620     }
621     }
622     }
623     else if (VTK_QUADRATIC_HEXAHEDRON==cellType)
624     {
625 dhawcroft 818 char tmpbuf2[20*20*2];
626 dhawcroft 796
627 dhawcroft 793 for (i = 0; i < numCells; i++)
628     {
629     if( elements->Id[i] >= elements->elementDistribution->vtxdist[myRank] && elements->Id[i] <= elements->elementDistribution->vtxdist[myRank+1]-1)
630     {
631     sprintf(tmpbuf2,"%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d\n",
632     nodesGlobal[elements->Nodes[INDEX2(0, i, NN)]],
633     nodesGlobal[elements->Nodes[INDEX2(1, i, NN)]],
634     nodesGlobal[elements->Nodes[INDEX2(2, i, NN)]],
635     nodesGlobal[elements->Nodes[INDEX2(3, i, NN)]],
636     nodesGlobal[elements->Nodes[INDEX2(4, i, NN)]],
637     nodesGlobal[elements->Nodes[INDEX2(5, i, NN)]],
638     nodesGlobal[elements->Nodes[INDEX2(6, i, NN)]],
639     nodesGlobal[elements->Nodes[INDEX2(7, i, NN)]],
640     nodesGlobal[elements->Nodes[INDEX2(8, i, NN)]],
641     nodesGlobal[elements->Nodes[INDEX2(9, i, NN)]],
642     nodesGlobal[elements->Nodes[INDEX2(10, i, NN)]],
643     nodesGlobal[elements->Nodes[INDEX2(11, i, NN)]],
644     nodesGlobal[elements->Nodes[INDEX2(16, i, NN)]],
645     nodesGlobal[elements->Nodes[INDEX2(17, i, NN)]],
646     nodesGlobal[elements->Nodes[INDEX2(18, i, NN)]],
647     nodesGlobal[elements->Nodes[INDEX2(19, i, NN)]],
648     nodesGlobal[elements->Nodes[INDEX2(12, i, NN)]],
649     nodesGlobal[elements->Nodes[INDEX2(13, i, NN)]],
650     nodesGlobal[elements->Nodes[INDEX2(14, i, NN)]],
651     nodesGlobal[elements->Nodes[INDEX2(15, i, NN)]]);
652 dhawcroft 818 __STRCAT(largebuf,tmpbuf2,tsz)
653 dhawcroft 793 elementCache.values[pos++]=i;
654     }
655     }
656     }
657     else if (numVTKNodesPerElement!=NN)
658     {
659 dhawcroft 796
660 dhawcroft 793 for (i = 0; i < numCells; i++)
661     {
662 dhawcroft 818 if (elements->Id[i] >= elements->elementDistribution->vtxdist[i] && elements->Id[i] <= elements->elementDistribution->vtxdist[i+1] - 1 )
663 dhawcroft 793 {
664 dhawcroft 818 for (j = 0; j < numVTKNodesPerElement; j++)
665     {
666     sprintf(tmpbuf,"%d ", nodesGlobal[elements->Nodes[INDEX2(elements->ReferenceElement->Type->geoNodes[j], i, NN)]]);
667     __STRCAT(largebuf,tmpbuf,tsz)
668     }
669     __STRCAT(largebuf,newline,tsz)
670     elementCache.values[pos++]=i;
671 dhawcroft 793 }
672     }
673     }
674     else
675 dhawcroft 796 {
676 dhawcroft 793 for(i = 0;i < numCells ; i++)
677     {
678     if( elements->Id[i] >= elements->elementDistribution->vtxdist[myRank] && elements->Id[i] <= elements->elementDistribution->vtxdist[myRank+1]-1)
679     {
680     for (j = 0; j < numVTKNodesPerElement; j++)
681     {
682     sprintf(tmpbuf,"%d ", nodesGlobal[ elements->Nodes[INDEX2(j, i, NN) ] ] );
683 dhawcroft 818 __STRCAT(largebuf,tmpbuf,tsz)
684 dhawcroft 793 }
685 dhawcroft 818 __STRCAT(largebuf,newline,tsz)
686 dhawcroft 793 elementCache.values[pos++]=i;
687     }
688     }
689 dhawcroft 796 }
690 dhawcroft 818
691 dhawcroft 796 elementCache.size = pos;
692 dhawcroft 793
693 dhawcroft 818 largebuf[tsz] = '\0';
694 dhawcroft 796 MPI_File_write_ordered(fh,largebuf,tsz, MPI_CHAR, &status);
695     MEMFREE(largebuf);
696 dhawcroft 793 MPIO_DEBUG(" Done Writing Connectivity ")
697 dhawcroft 794 MPIO_DEBUG(" Writing Offsets & Types... ")
698 dhawcroft 793
699     // Non-Collective
700     if( myRank == 0)
701     {
702     // write out the DataArray element for the offsets
703     char* tag1 = "</DataArray>\n<DataArray Name=\"offsets\" type=\"Int32\" format=\"ascii\">\n";
704     char* tag2 = "</DataArray>\n";
705     char *tag3 = "<DataArray Name=\"types\" type=\"UInt8\" format=\"ascii\">\n";
706     char *tag4 = "</DataArray>\n</Cells>\n";
707    
708     int n = numVTKNodesPerElement;
709    
710 dhawcroft 818 // allocate an upper bound on number of bytes needed
711 dhawcroft 793 int sz=0;
712     int lg = log10(numGlobalCells * n) + 1;
713 dhawcroft 818 sz += numGlobalCells*lg;
714     sz += numGlobalCells;
715 dhawcroft 796 tsz = 0;
716 dhawcroft 793
717     char* largebuf = MEMALLOC(sz + strlen(tag1) + strlen(tag2) + strlen(tag3) + strlen(tag4),char);
718     largebuf[0] ='\0';
719     char tmp[10];
720 dhawcroft 818 __STRCAT(largebuf,tag1,tsz)
721 dhawcroft 793 for (i=numVTKNodesPerElement; i<=numGlobalCells*numVTKNodesPerElement; i+=numVTKNodesPerElement)
722     {
723     sprintf(tmp,"%d\n", i);
724 dhawcroft 818 __STRCAT(largebuf,tmp,tsz)
725 dhawcroft 793 }
726 dhawcroft 818 __STRCAT(largebuf,tag2,tsz)
727     largebuf[tsz] = '\0';
728 dhawcroft 793 MPI_File_iwrite_shared(fh,largebuf, tsz,MPI_CHAR,&req);
729     MPI_Wait(&req,&status);
730 dhawcroft 818
731     // re-using buffer!
732 dhawcroft 793 largebuf[0] = '\0';
733     tsz = 0;
734 dhawcroft 818 __STRCAT(largebuf,tag3,tsz)
735 dhawcroft 793 for (i=0; i<numGlobalCells; i++)
736     {
737     sprintf(tmp, "%d\n", cellType);
738 dhawcroft 818 __STRCAT(largebuf,tmp,tsz)
739 dhawcroft 793 }
740 dhawcroft 818 __STRCAT(largebuf,tag4,tsz)
741     largebuf[tsz] = '\0';
742 dhawcroft 796 MPI_File_iwrite_shared(fh,largebuf,tsz,MPI_CHAR,&req);
743 dhawcroft 793 MPI_Wait(&req,&status);
744     MEMFREE(largebuf);
745     }
746    
747 dhawcroft 794 MPIO_DEBUG(" Done Writing Offsets & Types ")
748 dhawcroft 793
749 dhawcroft 818 // Write Cell data header Tags
750     if(myRank == 0)
751 dhawcroft 793 {
752 dhawcroft 818 MPIO_DEBUG(" Writing Cell Data ...")
753     if( write_celldata)
754 dhawcroft 793 {
755 dhawcroft 818 char tmpBuf[80];
756     char header[600];
757 dhawcroft 793 // mark the active data arrays
758     bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
759 dhawcroft 818 sprintf(tmpBuf, "<CellData");
760     strcat(header,tmpBuf);
761 dhawcroft 793 for (i_data =0 ;i_data<num_data;++i_data)
762     {
763 dhawcroft 818 if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data])
764 dhawcroft 793 {
765     // if the rank == 0: --> scalar data
766     // if the rank == 1: --> vector data
767     // if the rank == 2: --> tensor data
768    
769     switch(getDataPointRank(data_pp[i_data]))
770     {
771     case 0:
772     if (! set_scalar)
773     {
774     sprintf(tmpBuf," Scalars=\"%s\"",names_p[i_data]);
775     strcat(header,tmpBuf);
776     set_scalar=TRUE;
777     }
778     break;
779     case 1:
780     if (! set_vector)
781     {
782     sprintf(tmpBuf," Vectors=\"%s\"",names_p[i_data]);
783 dhawcroft 818 strcat(header,tmpBuf);
784 dhawcroft 793 set_vector=TRUE;
785     }
786     break;
787     case 2:
788     if (! set_tensor)
789     {
790     sprintf(tmpBuf," Tensors=\"%s\"",names_p[i_data]);
791 dhawcroft 818 strcat(header,tmpBuf);
792 dhawcroft 793 set_tensor=TRUE;
793     }
794     break;
795     default:
796     sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);
797     Finley_setError(VALUE_ERROR,error_msg);
798     return;
799     }
800     }
801     }
802     strcat(header, ">\n");
803     MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);
804     MPI_Wait(&req,&status);
805     }
806     }
807    
808 dhawcroft 818 // write actual data (collective)
809     if(write_celldata)
810 dhawcroft 793 {
811     for (i_data =0 ;i_data<num_data;++i_data)
812     {
813 dhawcroft 818 if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data])
814 dhawcroft 793 {
815     numPointsPerSample = elements->ReferenceElement->numQuadNodes;
816     rank = getDataPointRank(data_pp[i_data]);
817     nComp = getDataPointSize(data_pp[i_data]);
818     nCompReqd=1; // the number of components required by vtk
819     shape=0;
820     if (rank == 0)
821     {
822     nCompReqd = 1;
823     }
824     else if (rank == 1)
825     {
826     shape=getDataPointShape(data_pp[i_data], 0);
827     if (shape>3)
828     {
829     Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");
830     return;
831     }
832     nCompReqd = 3;
833     }
834     else
835     {
836     shape=getDataPointShape(data_pp[i_data], 0);
837     if (shape>3 || shape != getDataPointShape(data_pp[i_data], 1))
838     {
839     Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
840     return;
841     }
842     nCompReqd = 9;
843     }
844    
845     if( myRank == 0)
846     {
847     char header[250];
848 gross 850 sprintf(header,"<DataArray Name=\"%s\" type=\"Float64\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);
849 dhawcroft 793 MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);
850     MPI_Wait(&req,&status);
851     }
852    
853 dhawcroft 818 // Write the actual data
854     char tmpbuf[15];
855     char* largebuf = MEMALLOC(nCompReqd*numLocalCells*15 + numLocalCells*nCompReqd + nCompReqd + 15,char);
856 dhawcroft 793 largebuf[0] = '\0';
857     size_t tsz = 0;
858    
859 dhawcroft 818 double sampleAvg[nComp];
860    
861     for (k=0; k<elementCache.size; k++)
862 dhawcroft 793 {
863 dhawcroft 818 i = elementCache.values[k];
864 dhawcroft 793
865 dhawcroft 818 values = getSampleData(data_pp[i_data], i);
866     // averaging over the number of points in the sample
867     for (n=0; n<nComp; n++)
868 dhawcroft 793 {
869 gross 903 if (isExpanded(data_pp[i_data])) {
870     rtmp = 0.;
871     for (j=0; j<numPointsPerSample; j++) rtmp += values[INDEX2(n,j,nComp)];
872     sampleAvg[n] = rtmp/numPointsPerSample;
873     } else {
874     sampleAvg[n] = values[n];
875     }
876 dhawcroft 818 }
877     // if the number of required components is more than the number
878     // of actual components, pad with zeros
879    
880     // probably only need to get shape of first element
881     // write the data different ways for scalar, vector and tensor
882     if (nCompReqd == 1)
883     {
884     sprintf(tmpbuf, " %e", sampleAvg[0]);
885     __STRCAT(largebuf,tmpbuf,tsz)
886     }
887     else if (nCompReqd == 3)
888     {
889     // write out the data
890     for (m=0; m<shape; m++)
891 dhawcroft 793 {
892 dhawcroft 818 sprintf(tmpbuf, " %e", sampleAvg[m]);
893     __STRCAT(largebuf,tmpbuf,tsz)
894 dhawcroft 793 }
895 dhawcroft 818 for (m=0; m<nCompReqd-shape; m++)
896 dhawcroft 793 {
897 dhawcroft 818 __STRCAT(largebuf,zero,tsz)
898 dhawcroft 793 }
899     }
900 dhawcroft 818 else if (nCompReqd == 9)
901 dhawcroft 793 {
902 dhawcroft 818 // tensor data, so have a 3x3 matrix to output as a row
903     // of 9 data points
904     count = 0;
905     for (m=0; m<shape; m++)
906 dhawcroft 793 {
907 dhawcroft 818 for (n=0; n<shape; n++)
908 dhawcroft 793 {
909 dhawcroft 818 sprintf(tmpbuf, " %e", sampleAvg[count]);
910     __STRCAT(largebuf,tmpbuf,tsz)
911     count++;
912 dhawcroft 793 }
913 dhawcroft 818 for (n=0; n<3-shape; n++)
914 dhawcroft 793 {
915 dhawcroft 818 __STRCAT(largebuf,zero,tsz)
916 dhawcroft 793 }
917     }
918 dhawcroft 818 for (m=0; m<3-shape; m++)
919     for (n=0; n<3; n++)
920 dhawcroft 793 {
921 dhawcroft 818 __STRCAT(largebuf,zero,tsz)
922 dhawcroft 793 }
923     }
924 dhawcroft 818 __STRCAT(largebuf,newline,tsz)
925 dhawcroft 793 }
926 dhawcroft 818 largebuf[tsz] = '\0';
927 dhawcroft 793 MPI_File_write_ordered(fh,largebuf,tsz,MPI_CHAR,&status);
928     MEMFREE(largebuf);
929     if( myRank == 0)
930     {
931     char *tag = "</DataArray>\n";
932     MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);
933     MPI_Wait(&req,&status);
934     }
935 dhawcroft 818
936 dhawcroft 793 }
937     }
938 dhawcroft 818 // closing celldata tag
939 dhawcroft 793 if(myRank == 0)
940     {
941 dhawcroft 818 char* tag = "</CellData>\n";
942 dhawcroft 793 MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);
943     MPI_Wait(&req,&status);
944     }
945 dhawcroft 818
946     MPIO_DEBUG(" Done Writing Cell Data ")
947 dhawcroft 793 }
948    
949 dhawcroft 818
950     // Write Point Data Header Tags
951     if( myRank == 0)
952 dhawcroft 793 {
953 dhawcroft 818 char header[600];
954     char tmpBuf[50];
955    
956     if (write_pointdata)
957 dhawcroft 793 {
958 dhawcroft 818 MPIO_DEBUG(" Writing Pointdata... ")
959 dhawcroft 793 // mark the active data arrays
960     bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
961 dhawcroft 818 sprintf(header, "<PointData");
962 dhawcroft 793 for (i_data =0 ;i_data<num_data;++i_data)
963     {
964 dhawcroft 818 if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])
965 dhawcroft 793 {
966     // if the rank == 0: --> scalar data
967     // if the rank == 1: --> vector data
968     // if the rank == 2: --> tensor data
969    
970     switch(getDataPointRank(data_pp[i_data]))
971     {
972     case 0:
973     if (! set_scalar)
974     {
975     sprintf(tmpBuf," Scalars=\"%s\"",names_p[i_data]);
976     strcat(header,tmpBuf);
977     set_scalar=TRUE;
978     }
979     break;
980     case 1:
981     if (! set_vector)
982     {
983     sprintf(tmpBuf," Vectors=\"%s\"",names_p[i_data]);
984 dhawcroft 818 strcat(header,tmpBuf);
985 dhawcroft 793 set_vector=TRUE;
986     }
987     break;
988     case 2:
989     if (! set_tensor)
990     {
991     sprintf(tmpBuf," Tensors=\"%s\"",names_p[i_data]);
992 dhawcroft 818 strcat(header,tmpBuf);
993 dhawcroft 793 set_tensor=TRUE;
994     }
995     break;
996     default:
997     sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);
998     Finley_setError(VALUE_ERROR,error_msg);
999     return;
1000     }
1001     }
1002     }
1003     strcat(header, ">\n");
1004     MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);
1005     MPI_Wait(&req,&status);
1006     }
1007     }
1008    
1009 dhawcroft 818 // write actual data
1010     if(write_pointdata)
1011 dhawcroft 793 {
1012     for (i_data =0 ;i_data<num_data;++i_data)
1013     {
1014 dhawcroft 818 if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])
1015 dhawcroft 793 {
1016     numPointsPerSample = elements->ReferenceElement->numQuadNodes;
1017     rank = getDataPointRank(data_pp[i_data]);
1018     nComp = getDataPointSize(data_pp[i_data]);
1019     nCompReqd=1; // the number of components required by vtk
1020     shape=0;
1021     if (rank == 0)
1022     {
1023     nCompReqd = 1;
1024     }
1025     else if (rank == 1)
1026     {
1027     shape=getDataPointShape(data_pp[i_data], 0);
1028     if (shape>3)
1029     {
1030     Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");
1031     return;
1032     }
1033     nCompReqd = 3;
1034     }
1035     else
1036     {
1037     shape=getDataPointShape(data_pp[i_data], 0);
1038     if (shape>3 || shape != getDataPointShape(data_pp[i_data], 1))
1039     {
1040     Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
1041     return;
1042     }
1043     nCompReqd = 9;
1044     }
1045    
1046     if( myRank == 0)
1047     {
1048     char header[250];
1049 gross 850 sprintf(header, "<DataArray Name=\"%s\" type=\"Float64\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);
1050 dhawcroft 793 MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);
1051     MPI_Wait(&req,&status);
1052     }
1053 dhawcroft 818 // write out the data
1054     // if the number of required components is more than the number
1055     // of actual components, pad with zeros
1056 dhawcroft 793
1057 dhawcroft 818 char tmpbuf[15];
1058     char* largebuf = MEMALLOC(nCompReqd*numLocalNodes*15 + numLocal*nCompReqd + nCompReqd + 15,char);
1059 dhawcroft 793 largebuf[0] = '\0';
1060 dhawcroft 818 bool_t do_write=TRUE;
1061 dhawcroft 793 size_t tsz = 0;
1062    
1063 dhawcroft 818 for(k=0;k < nodeCache.size;k++)
1064 dhawcroft 793 {
1065 dhawcroft 818 i = nodeCache.values[k];
1066 dhawcroft 793
1067 dhawcroft 818 if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1068 dhawcroft 793 {
1069 dhawcroft 818 if (mesh_p->Nodes->toReduced[i]>=0)
1070     {
1071     switch(getFunctionSpaceType(data_pp[i_data]))
1072     {
1073     case FINLEY_DEGREES_OF_FREEDOM:
1074     values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);
1075     break;
1076     case FINLEY_REDUCED_DEGREES_OF_FREEDOM:
1077     values = getSampleData(data_pp[i_data],mesh_p->Nodes->reducedDegreeOfFreedom[i]);
1078     break;
1079     case FINLEY_NODES:
1080     values = getSampleData(data_pp[i_data],i);
1081     break;
1082     }
1083     do_write=TRUE;
1084     }
1085     else
1086     {
1087     do_write=FALSE;
1088     }
1089 dhawcroft 793 }
1090 dhawcroft 818 else
1091 dhawcroft 793 {
1092 dhawcroft 818 do_write=TRUE;
1093     switch(getFunctionSpaceType(data_pp[i_data]))
1094     {
1095     case FINLEY_DEGREES_OF_FREEDOM:
1096     values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);
1097     break;
1098     case FINLEY_NODES:
1099     values = getSampleData(data_pp[i_data],i);
1100     break;
1101     }
1102 dhawcroft 793 }
1103 dhawcroft 818 // write the data different ways for scalar, vector and tensor
1104     if (do_write)
1105 dhawcroft 793 {
1106 dhawcroft 818 if (nCompReqd == 1)
1107 dhawcroft 793 {
1108 dhawcroft 818 sprintf(tmpbuf," %e", values[0]);
1109     __STRCAT(largebuf,tmpbuf,tsz)
1110 dhawcroft 793 }
1111 dhawcroft 818 else if (nCompReqd == 3)
1112 dhawcroft 793 {
1113 dhawcroft 818 for (m=0; m<shape; m++)
1114     {
1115 dhawcroft 793
1116 dhawcroft 818 sprintf(tmpbuf," %e",values[m]);
1117     __STRCAT(largebuf,tmpbuf,tsz)
1118     }
1119     for (m=0; m<nCompReqd-shape; m++)
1120     {
1121     __STRCAT(largebuf,zero,tsz)
1122     }
1123 dhawcroft 793 }
1124 dhawcroft 818 else if (nCompReqd == 9)
1125 dhawcroft 793 {
1126 dhawcroft 818 // tensor data, so have a 3x3 matrix to output as a row
1127     // of 9 data points
1128     count = 0;
1129     for (m=0; m<shape; m++)
1130 dhawcroft 793 {
1131 dhawcroft 818 for (n=0; n<shape; n++)
1132     {
1133     sprintf(tmpbuf," %e", values[count]);
1134     __STRCAT(largebuf,tmpbuf,tsz)
1135     count++;
1136     }
1137     for (n=0; n<3-shape; n++)
1138     {
1139     __STRCAT(largebuf,zero,tsz)
1140     }
1141 dhawcroft 793 }
1142 dhawcroft 818 for (m=0; m<3-shape; m++)
1143 dhawcroft 793 {
1144 dhawcroft 818 for (n=0; n<3; n++)
1145     {
1146     __STRCAT(largebuf,zero,tsz)
1147     }
1148 dhawcroft 793 }
1149     }
1150 dhawcroft 818 __STRCAT(largebuf,newline,tsz)
1151 dhawcroft 793 }
1152 dhawcroft 818
1153 dhawcroft 793 }
1154 dhawcroft 818 // Write out local data
1155    
1156     largebuf[tsz] = '\0';
1157 dhawcroft 793 MPI_File_write_ordered(fh,largebuf,tsz,MPI_CHAR,&status);
1158     MEMFREE(largebuf);
1159     if( myRank == 0)
1160     {
1161     char *tag = "</DataArray>\n";
1162     MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);
1163     MPI_Wait(&req,&status);
1164     }
1165     }
1166     }
1167 dhawcroft 818 // Finish off with closing tag
1168 dhawcroft 793 if(myRank == 0)
1169     {
1170 dhawcroft 818 char* tag = "</PointData>\n";
1171 dhawcroft 793 MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);
1172     MPI_Wait(&req,&status);
1173     }
1174 dhawcroft 818 MPIO_DEBUG(" Done Writing Pointdata ")
1175 dhawcroft 793 }
1176 dhawcroft 818 // end write_pointdata
1177 dhawcroft 793
1178 dhawcroft 818 // tag and bag...
1179 dhawcroft 793 if (myRank == 0)
1180     {
1181     char *footer = "</Piece>\n</UnstructuredGrid>\n</VTKFile>";
1182     MPI_File_iwrite_shared(fh,footer,strlen(footer),MPI_CHAR,&req);
1183     MPI_Wait(&req,&status);
1184     }
1185    
1186     MEMFREE(nodesGlobal);
1187     MEMFREE(nodeCache.values);
1188     MEMFREE(elementCache.values);
1189     #ifdef MPIO_HINTS
1190     MPI_Info_free(&infoHints);
1191 dhawcroft 796 #undef MPIO_HINTS
1192 dhawcroft 793 #endif
1193     MPI_File_close(&fh);
1194     MPIO_DEBUG(" ***** Exit saveVTK ***** ")
1195 dhawcroft 796 #undef __STRCAT
1196 dhawcroft 793 }
1197    
1198     #undef MPIO_DEBUG
1199     #else
1200    
1201    
1202    
1203    
1204     void Finley_Mesh_saveVTK(const char * filename_p, Finley_Mesh *mesh_p, const dim_t num_data,char* *names_p, escriptDataC* *data_pp)
1205     {
1206 gross 1028 #define NCOMP_MAX 9
1207 dhawcroft 793 char error_msg[LenErrorMsg_MAX];
1208 gross 1028 double sampleAvg[NCOMP_MAX];
1209 jgs 110 /* if there is no mesh we just return */
1210 jgs 153
1211 gross 687 int i, j, k, numVTKNodesPerElement,i_data,m, count, n, rank,shape, numPoints, cellType, numCells,
1212 gross 1028 nDim, numPointsPerSample, nComp, nCompReqd, NN;
1213 jgs 150 index_t j2;
1214 gross 1028 int nodetype=FINLEY_DEGREES_OF_FREEDOM;
1215     int elementtype=FINLEY_UNKNOWN;
1216 jgs 113 double* values, rtmp;
1217 gross 687 char elemTypeStr[32];
1218 gross 1028 FILE * fileHandle_p = NULL;
1219     bool_t do_write, *isCellCentered=NULL,write_celldata=FALSE,write_pointdata=FALSE;
1220     Finley_ElementFile* elements=NULL;
1221     ElementTypeId TypeId;
1222 dhawcroft 793
1223 jgs 153 /* open the file and check handle */
1224 gross 1028 if (mesh_p==NULL) return;
1225 jgs 153
1226 gross 1028 fileHandle_p = fopen(filename_p, "w");
1227 dhawcroft 793 if (fileHandle_p==NULL)
1228     {
1229 jgs 153 sprintf(error_msg, "saveVTK: File %s could not be opened for writing.", filename_p);
1230     Finley_setError(IO_ERROR,error_msg);
1231     return;
1232 jgs 150 }
1233 jgs 153 /* find the mesh type to be written */
1234 gross 1028 isCellCentered=TMPMEMALLOC(num_data,bool_t);
1235    
1236    
1237     if (Finley_checkPtr(isCellCentered)) {
1238     fclose(fileHandle_p);
1239     return;
1240     }
1241 dhawcroft 793 for (i_data=0;i_data<num_data;++i_data)
1242     {
1243     if (! isEmpty(data_pp[i_data]))
1244     {
1245     switch(getFunctionSpaceType(data_pp[i_data]))
1246     {
1247     case FINLEY_DEGREES_OF_FREEDOM:
1248     nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1249     if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)
1250     {
1251     elementtype=FINLEY_ELEMENTS;
1252 jgs 153 }
1253 dhawcroft 793 else
1254     {
1255     Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1256     fclose(fileHandle_p);
1257 gross 1028 TMPMEMFREE(isCellCentered);
1258 dhawcroft 793 return;
1259 jgs 153 }
1260 dhawcroft 793 isCellCentered[i_data]=FALSE;
1261     break;
1262     case FINLEY_REDUCED_DEGREES_OF_FREEDOM:
1263     nodetype = FINLEY_REDUCED_DEGREES_OF_FREEDOM;
1264     if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)
1265     {
1266     elementtype=FINLEY_ELEMENTS;
1267     }
1268     else
1269     {
1270     Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1271     fclose(fileHandle_p);
1272 gross 1028 TMPMEMFREE(isCellCentered);
1273 dhawcroft 793 return;
1274     }
1275     isCellCentered[i_data]=FALSE;
1276     break;
1277     case FINLEY_NODES:
1278     nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1279     if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)
1280     {
1281     elementtype=FINLEY_ELEMENTS;
1282     }
1283     else
1284     {
1285     Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1286     fclose(fileHandle_p);
1287 gross 1028 TMPMEMFREE(isCellCentered);
1288 dhawcroft 793 return;
1289     }
1290     isCellCentered[i_data]=FALSE;
1291     break;
1292     case FINLEY_ELEMENTS:
1293     nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1294     if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)
1295     {
1296     elementtype=FINLEY_ELEMENTS;
1297     }
1298     else
1299     {
1300     Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1301     fclose(fileHandle_p);
1302     return;
1303 gross 1028 TMPMEMFREE(isCellCentered);
1304 dhawcroft 793 }
1305     isCellCentered[i_data]=TRUE;
1306     break;
1307     case FINLEY_FACE_ELEMENTS:
1308     nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1309     if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_FACE_ELEMENTS)
1310     {
1311     elementtype=FINLEY_FACE_ELEMENTS;
1312     }
1313     else
1314     {
1315     Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1316     fclose(fileHandle_p);
1317 gross 1028 TMPMEMFREE(isCellCentered);
1318 dhawcroft 793 return;
1319     }
1320     isCellCentered[i_data]=TRUE;
1321     break;
1322     case FINLEY_POINTS:
1323     nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1324     if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_POINTS)
1325     {
1326     elementtype=FINLEY_POINTS;
1327     }
1328     else
1329     {
1330     Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1331     fclose(fileHandle_p);
1332 gross 1028 TMPMEMFREE(isCellCentered);
1333 dhawcroft 793 return;
1334     }
1335     isCellCentered[i_data]=TRUE;
1336     break;
1337     case FINLEY_CONTACT_ELEMENTS_1:
1338     nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1339     if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1)
1340     {
1341     elementtype=FINLEY_CONTACT_ELEMENTS_1;
1342     }
1343     else
1344     {
1345     Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1346     fclose(fileHandle_p);
1347 gross 1028 TMPMEMFREE(isCellCentered);
1348 dhawcroft 793 return;
1349     }
1350     isCellCentered[i_data]=TRUE;
1351     break;
1352     case FINLEY_CONTACT_ELEMENTS_2:
1353     nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1354     if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1)
1355     {
1356     elementtype=FINLEY_CONTACT_ELEMENTS_1;
1357     }
1358     else
1359     {
1360     Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1361     fclose(fileHandle_p);
1362 gross 1028 TMPMEMFREE(isCellCentered);
1363 dhawcroft 793 return;
1364     }
1365     isCellCentered[i_data]=TRUE;
1366     break;
1367     default:
1368     sprintf(error_msg,"saveVTK: Finley does not know anything about function space type %d",getFunctionSpaceType(data_pp[i_data]));
1369     Finley_setError(TYPE_ERROR,error_msg);
1370     fclose(fileHandle_p);
1371 gross 1028 TMPMEMFREE(isCellCentered);
1372 dhawcroft 793 return;
1373     }
1374     if (isCellCentered[i_data])
1375     {
1376     write_celldata=TRUE;
1377     }
1378     else
1379     {
1380     write_pointdata=TRUE;
1381     }
1382     }
1383 jgs 153 }
1384     /* select nomber of points and the mesh component */
1385 gross 687 numPoints = mesh_p->Nodes->numNodes;
1386 dhawcroft 793 if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1387     {
1388     numPoints = mesh_p->Nodes->reducedNumNodes;
1389 jgs 153 }
1390 dhawcroft 793 else
1391     {
1392     numPoints = mesh_p->Nodes->numNodes;
1393     }
1394 jgs 153 if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;
1395 dhawcroft 793 switch(elementtype)
1396     {
1397     case FINLEY_ELEMENTS:
1398     elements=mesh_p->Elements;
1399     break;
1400     case FINLEY_FACE_ELEMENTS:
1401     elements=mesh_p->FaceElements;
1402     break;
1403     case FINLEY_POINTS:
1404     elements=mesh_p->Points;
1405     break;
1406     case FINLEY_CONTACT_ELEMENTS_1:
1407     elements=mesh_p->ContactElements;
1408     break;
1409 jgs 113 }
1410 dhawcroft 793 if (elements==NULL)
1411     {
1412     Finley_setError(SYSTEM_ERROR,"saveVTK: undefined element file");
1413     fclose(fileHandle_p);
1414 gross 1028 TMPMEMFREE(isCellCentered);
1415 dhawcroft 793 return;
1416 jgs 113 }
1417 jgs 153 /* map finley element type to VTK element type */
1418 dhawcroft 793 numCells = elements->numElements;
1419     if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1420     {
1421     TypeId = elements->LinearReferenceElement->Type->TypeId;
1422 jgs 113 }
1423 dhawcroft 793 else
1424     {
1425     TypeId = elements->ReferenceElement->Type->TypeId;
1426     }
1427 jgs 110
1428 dhawcroft 793 switch(TypeId)
1429     {
1430     case Point1:
1431     case Line2Face:
1432     case Line3Face:
1433     case Point1_Contact:
1434     case Line2Face_Contact:
1435     case Line3Face_Contact:
1436     cellType = VTK_VERTEX;
1437     numVTKNodesPerElement = 1;
1438     strcpy(elemTypeStr, "VTK_VERTEX");
1439     break;
1440 jgs 153
1441 dhawcroft 793 case Line2:
1442     case Tri3Face:
1443     case Rec4Face:
1444     case Line2_Contact:
1445     case Tri3_Contact:
1446     case Tri3Face_Contact:
1447     case Rec4Face_Contact:
1448     cellType = VTK_LINE;
1449     numVTKNodesPerElement = 2;
1450     strcpy(elemTypeStr, "VTK_LINE");
1451     break;
1452 jgs 153
1453 dhawcroft 793 case Tri3:
1454     case Tet4Face:
1455     case Tet4Face_Contact:
1456     cellType = VTK_TRIANGLE;
1457     numVTKNodesPerElement = 3;
1458     strcpy(elemTypeStr, "VTK_TRIANGLE");
1459     break;
1460 jgs 153
1461 dhawcroft 793 case Rec4:
1462     case Hex8Face:
1463     case Rec4_Contact:
1464     case Hex8Face_Contact:
1465     cellType = VTK_QUAD;
1466     numVTKNodesPerElement = 4;
1467     strcpy(elemTypeStr, "VTK_QUAD");
1468     break;
1469 jgs 153
1470 dhawcroft 793 case Tet4:
1471     cellType = VTK_TETRA;
1472     numVTKNodesPerElement = 4;
1473     strcpy(elemTypeStr, "VTK_TETRA");
1474     break;
1475 jgs 153
1476 dhawcroft 793 case Hex8:
1477     cellType = VTK_HEXAHEDRON;
1478     numVTKNodesPerElement = 8;
1479     strcpy(elemTypeStr, "VTK_HEXAHEDRON");
1480     break;
1481 jgs 153
1482 dhawcroft 793 case Line3:
1483     case Tri6Face:
1484     case Rec8Face:
1485     case Line3_Contact:
1486     case Tri6Face_Contact:
1487     case Rec8Face_Contact:
1488     cellType = VTK_QUADRATIC_EDGE;
1489     numVTKNodesPerElement = 3;
1490     strcpy(elemTypeStr, "VTK_QUADRATIC_EDGE");
1491     break;
1492 jgs 153
1493 dhawcroft 793 case Tri6:
1494     case Tet10Face:
1495     case Tri6_Contact:
1496     case Tet10Face_Contact:
1497     cellType = VTK_QUADRATIC_TRIANGLE;
1498     numVTKNodesPerElement = 6;
1499     strcpy(elemTypeStr, "VTK_QUADRATIC_TRIANGLE");
1500     break;
1501 jgs 153
1502 dhawcroft 793 case Rec8:
1503     case Hex20Face:
1504     case Rec8_Contact:
1505     case Hex20Face_Contact:
1506     cellType = VTK_QUADRATIC_QUAD;
1507     numVTKNodesPerElement = 8;
1508     strcpy(elemTypeStr, "VTK_QUADRATIC_QUAD");
1509     break;
1510 jgs 153
1511 dhawcroft 793 case Tet10:
1512     cellType = VTK_QUADRATIC_TETRA;
1513     numVTKNodesPerElement = 10;
1514     strcpy(elemTypeStr, "VTK_QUADRATIC_TETRA");
1515     break;
1516 jgs 153
1517 dhawcroft 793 case Hex20:
1518     cellType = VTK_QUADRATIC_HEXAHEDRON;
1519     numVTKNodesPerElement = 20;
1520     strcpy(elemTypeStr, "VTK_QUADRATIC_HEXAHEDRON");
1521     break;
1522 jgs 153
1523 dhawcroft 793 default:
1524     sprintf(error_msg, "saveVTK: Element type %s is not supported by VTK",elements->ReferenceElement->Type->Name);
1525     Finley_setError(VALUE_ERROR,error_msg);
1526     fclose(fileHandle_p);
1527 gross 1028 TMPMEMFREE(isCellCentered);
1528 dhawcroft 793 return;
1529     }
1530 jgs 153 /* xml header */
1531     fprintf(fileHandle_p, "<?xml version=\"1.0\"?>\n");
1532     fprintf(fileHandle_p, "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\">\n");
1533 jgs 110
1534 jgs 153 /* finley uses an unstructured mesh, so UnstructuredGrid *should* work */
1535     fprintf(fileHandle_p, "<UnstructuredGrid>\n");
1536    
1537     /* is there only one "piece" to the data?? */
1538     fprintf(fileHandle_p, "<Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n",numPoints, numCells);
1539     /* now for the points; equivalent to positions section in saveDX() */
1540     /* "The points element explicitly defines coordinates for each point
1541 dhawcroft 793 * individually. It contains one DataArray element describing an array
1542     * with three components per value, each specifying the coordinates of one
1543     * point" - from Vtk User's Guide
1544     */
1545 jgs 153 fprintf(fileHandle_p, "<Points>\n");
1546 dhawcroft 793 /*
1547     * the reason for this if statement is explained in the long comment below
1548     */
1549 gross 687 nDim = mesh_p->Nodes->numDim;
1550 gross 850 fprintf(fileHandle_p, "<DataArray NumberOfComponents=\"%d\" type=\"Float64\" format=\"ascii\">\n",MAX(3,nDim));
1551 jgs 153 /* vtk/mayavi doesn't like 2D data, it likes 3D data with a degenerate
1552 dhawcroft 793 * third dimension to handle 2D data (like a sheet of paper). So, if
1553     * nDim is 2, we have to append zeros to the array to get this third
1554     * dimension, and keep the visualisers happy.
1555     * Indeed, if nDim is less than 3, must pad all empty dimensions, so
1556     * that the total number of dims is 3.
1557 jgs 153 */
1558 dhawcroft 793 if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1559     {
1560     for (i = 0; i < mesh_p->Nodes->numNodes; i++)
1561     {
1562     if (mesh_p->Nodes->toReduced[i]>=0)
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     else
1571     {
1572     for (i = 0; i < mesh_p->Nodes->numNodes; i++)
1573     {
1574    
1575     for (j = 0; j < nDim; j++) fprintf(fileHandle_p, " %e",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);
1576     for (k=0; k<3-nDim; k++) fprintf(fileHandle_p, " %e",0.);
1577     fprintf(fileHandle_p, "\n");
1578     }
1579    
1580     }
1581 jgs 153 fprintf(fileHandle_p, "</DataArray>\n");
1582     fprintf(fileHandle_p, "</Points>\n");
1583 jgs 113
1584 jgs 153 /* write out the DataArray element for the connectivity */
1585 jgs 113
1586 gross 1028 NN = elements->ReferenceElement->Type->numNodes;
1587 jgs 153 fprintf(fileHandle_p, "<Cells>\n");
1588     fprintf(fileHandle_p, "<DataArray Name=\"connectivity\" type=\"Int32\" format=\"ascii\">\n");
1589 jgs 113
1590 dhawcroft 793 if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1591     {
1592     for (i = 0; i < numCells; i++)
1593     {
1594     for (j = 0; j < numVTKNodesPerElement; j++)
1595     fprintf(fileHandle_p,"%d ",mesh_p->Nodes->toReduced[elements->Nodes[INDEX2(elements->ReferenceElement->Type->linearNodes[j], i, NN)]]);
1596     fprintf(fileHandle_p, "\n");
1597     }
1598     }
1599     else if (VTK_QUADRATIC_HEXAHEDRON==cellType)
1600     {
1601     for (i = 0; i < numCells; i++)
1602     {
1603     fprintf(fileHandle_p,"%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d\n",
1604     elements->Nodes[INDEX2(0, i, NN)],
1605     elements->Nodes[INDEX2(1, i, NN)],
1606     elements->Nodes[INDEX2(2, i, NN)],
1607     elements->Nodes[INDEX2(3, i, NN)],
1608     elements->Nodes[INDEX2(4, i, NN)],
1609     elements->Nodes[INDEX2(5, i, NN)],
1610     elements->Nodes[INDEX2(6, i, NN)],
1611     elements->Nodes[INDEX2(7, i, NN)],
1612     elements->Nodes[INDEX2(8, i, NN)],
1613     elements->Nodes[INDEX2(9, i, NN)],
1614     elements->Nodes[INDEX2(10, i, NN)],
1615     elements->Nodes[INDEX2(11, i, NN)],
1616     elements->Nodes[INDEX2(16, i, NN)],
1617     elements->Nodes[INDEX2(17, i, NN)],
1618     elements->Nodes[INDEX2(18, i, NN)],
1619     elements->Nodes[INDEX2(19, i, NN)],
1620     elements->Nodes[INDEX2(12, i, NN)],
1621     elements->Nodes[INDEX2(13, i, NN)],
1622     elements->Nodes[INDEX2(14, i, NN)],
1623     elements->Nodes[INDEX2(15, i, NN)]);
1624     }
1625     }
1626     else if (numVTKNodesPerElement!=NN)
1627     {
1628     for (i = 0; i < numCells; i++)
1629     {
1630     for (j = 0; j < numVTKNodesPerElement; j++) fprintf(fileHandle_p,"%d ", elements->Nodes[INDEX2(elements->ReferenceElement->Type->geoNodes[j], i, NN)]);
1631     fprintf(fileHandle_p, "\n");
1632     }
1633     }
1634     else
1635     {
1636     for (i = 0; i < numCells; i++)
1637     {
1638     for (j = 0; j < numVTKNodesPerElement; j++) fprintf(fileHandle_p,"%d ", elements->Nodes[INDEX2(j, i, NN)]);
1639     fprintf(fileHandle_p, "\n");
1640     }
1641     }
1642 jgs 153 fprintf(fileHandle_p, "</DataArray>\n");
1643 jgs 110
1644 jgs 153 /* write out the DataArray element for the offsets */
1645     fprintf(fileHandle_p, "<DataArray Name=\"offsets\" type=\"Int32\" format=\"ascii\">\n");
1646     for (i=numVTKNodesPerElement; i<=numCells*numVTKNodesPerElement; i+=numVTKNodesPerElement) fprintf(fileHandle_p, "%d\n", i);
1647     fprintf(fileHandle_p, "</DataArray>\n");
1648    
1649     /* write out the DataArray element for the types */
1650     fprintf(fileHandle_p, "<DataArray Name=\"types\" type=\"UInt8\" format=\"ascii\">\n");
1651     for (i=0; i<numCells; i++) fprintf(fileHandle_p, "%d\n", cellType);
1652     fprintf(fileHandle_p, "</DataArray>\n");
1653    
1654     /* finish off the <Cells> element */
1655     fprintf(fileHandle_p, "</Cells>\n");
1656    
1657     /* cell data */
1658 dhawcroft 793 if (write_celldata)
1659     {
1660     /* mark the active data arrays */
1661     bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
1662     fprintf(fileHandle_p, "<CellData");
1663     for (i_data =0 ;i_data<num_data;++i_data)
1664     {
1665     if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data])
1666     {
1667     /* if the rank == 0: --> scalar data
1668     * if the rank == 1: --> vector data
1669     * if the rank == 2: --> tensor data
1670     */
1671     switch(getDataPointRank(data_pp[i_data]))
1672     {
1673     case 0:
1674     if (! set_scalar)
1675     {
1676     fprintf(fileHandle_p," Scalars=\"%s\"",names_p[i_data]);
1677     set_scalar=TRUE;
1678     }
1679     break;
1680     case 1:
1681     if (! set_vector)
1682     {
1683     fprintf(fileHandle_p," Vectors=\"%s\"",names_p[i_data]);
1684     set_vector=TRUE;
1685     }
1686     break;
1687     case 2:
1688     if (! set_tensor)
1689     {
1690     fprintf(fileHandle_p," Tensors=\"%s\"",names_p[i_data]);
1691     set_tensor=TRUE;
1692     }
1693     break;
1694     default:
1695     sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);
1696     Finley_setError(VALUE_ERROR,error_msg);
1697     fclose(fileHandle_p);
1698 gross 1028 TMPMEMFREE(isCellCentered);
1699 dhawcroft 793 return;
1700     }
1701     }
1702     }
1703     fprintf(fileHandle_p, ">\n");
1704     /* write the arrays */
1705     for (i_data =0 ;i_data<num_data;++i_data)
1706     {
1707     if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data])
1708     {
1709     numPointsPerSample = elements->ReferenceElement->numQuadNodes;
1710     rank = getDataPointRank(data_pp[i_data]);
1711     nComp = getDataPointSize(data_pp[i_data]);
1712     nCompReqd=1; /* the number of components required by vtk */
1713     shape=0;
1714     if (rank == 0)
1715     {
1716     nCompReqd = 1;
1717     }
1718     else if (rank == 1)
1719     {
1720     shape=getDataPointShape(data_pp[i_data], 0);
1721     if (shape>3)
1722     {
1723     Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");
1724     fclose(fileHandle_p);
1725 gross 1028 TMPMEMFREE(isCellCentered);
1726 dhawcroft 793 return;
1727     }
1728     nCompReqd = 3;
1729     }
1730     else
1731     {
1732     shape=getDataPointShape(data_pp[i_data], 0);
1733     if (shape>3 || shape != getDataPointShape(data_pp[i_data], 1))
1734     {
1735     Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
1736     fclose(fileHandle_p);
1737 gross 1028 TMPMEMFREE(isCellCentered);
1738 dhawcroft 793 return;
1739     }
1740     nCompReqd = 9;
1741     }
1742 gross 850 fprintf(fileHandle_p, "<DataArray Name=\"%s\" type=\"Float64\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);
1743 dhawcroft 793
1744     for (i=0; i<numCells; i++)
1745     {
1746     values = getSampleData(data_pp[i_data], i);
1747     /* averaging over the number of points in the sample */
1748 gross 1028 for (k=0; k<MIN(nComp,NCOMP_MAX); k++)
1749 dhawcroft 793 {
1750 gross 903 if (isExpanded(data_pp[i_data])) {
1751     rtmp = 0.;
1752     for (j=0; j<numPointsPerSample; j++) rtmp += values[INDEX2(k,j,nComp)];
1753     sampleAvg[k] = rtmp/numPointsPerSample;
1754     } else {
1755     sampleAvg[k] = values[k];
1756     }
1757    
1758 dhawcroft 793 }
1759     /* if the number of required components is more than the number
1760     * of actual components, pad with zeros
1761     */
1762     /* probably only need to get shape of first element */
1763     /* write the data different ways for scalar, vector and tensor */
1764     if (nCompReqd == 1)
1765     {
1766     fprintf(fileHandle_p, " %e", sampleAvg[0]);
1767     }
1768     else if (nCompReqd == 3)
1769     {
1770     /* write out the data */
1771     for (m=0; m<shape; m++) fprintf(fileHandle_p, " %e", sampleAvg[m]);
1772     for (m=0; m<nCompReqd-shape; m++) fprintf(fileHandle_p, " %e", 0.);
1773     }
1774     else if (nCompReqd == 9)
1775     {
1776     /* tensor data, so have a 3x3 matrix to output as a row
1777     * of 9 data points */
1778     count = 0;
1779     for (m=0; m<shape; m++)
1780     {
1781     for (n=0; n<shape; n++)
1782     {
1783     fprintf(fileHandle_p, " %e", sampleAvg[count]);
1784     count++;
1785     }
1786     for (n=0; n<3-shape; n++) fprintf(fileHandle_p, " %e", 0.);
1787 jgs 150 }
1788 dhawcroft 793 for (m=0; m<3-shape; m++)
1789     for (n=0; n<3; n++) fprintf(fileHandle_p, " %e", 0.);
1790     }
1791     fprintf(fileHandle_p, "\n");
1792     }
1793     fprintf(fileHandle_p, "</DataArray>\n");
1794     }
1795     }
1796     fprintf(fileHandle_p, "</CellData>\n");
1797 jgs 110 }
1798 jgs 153 /* point data */
1799 dhawcroft 793 if (write_pointdata)
1800     {
1801     /* mark the active data arrays */
1802     bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
1803     fprintf(fileHandle_p, "<PointData");
1804     for (i_data =0 ;i_data<num_data;++i_data)
1805     {
1806     if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])
1807     {
1808     /* if the rank == 0: --> scalar data
1809     * if the rank == 1: --> vector data
1810     * if the rank == 2: --> tensor data
1811     */
1812     switch(getDataPointRank(data_pp[i_data]))
1813     {
1814     case 0:
1815     if (! set_scalar)
1816     {
1817     fprintf(fileHandle_p," Scalars=\"%s\"",names_p[i_data]);
1818     set_scalar=TRUE;
1819     }
1820     break;
1821     case 1:
1822     if (! set_vector)
1823     {
1824     fprintf(fileHandle_p," Vectors=\"%s\"",names_p[i_data]);
1825     set_vector=TRUE;
1826     }
1827     break;
1828     case 2:
1829     if (! set_tensor)
1830     {
1831     fprintf(fileHandle_p," Tensors=\"%s\"",names_p[i_data]);
1832     set_tensor=TRUE;
1833     }
1834     break;
1835     default:
1836     sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);
1837     Finley_setError(VALUE_ERROR,error_msg);
1838     fclose(fileHandle_p);
1839 gross 1028 TMPMEMFREE(isCellCentered);
1840 dhawcroft 793 return;
1841     }
1842     }
1843     }
1844     fprintf(fileHandle_p, ">\n");
1845     /* write the arrays */
1846     for (i_data =0 ;i_data<num_data;++i_data)
1847     {
1848     if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])
1849     {
1850     numPointsPerSample = elements->ReferenceElement->numQuadNodes;
1851     rank = getDataPointRank(data_pp[i_data]);
1852     nComp = getDataPointSize(data_pp[i_data]);
1853     nCompReqd=1; /* the number of components required by vtk */
1854     shape=0;
1855     if (rank == 0)
1856     {
1857     nCompReqd = 1;
1858     }
1859     else if (rank == 1)
1860     {
1861     shape=getDataPointShape(data_pp[i_data], 0);
1862     if (shape>3)
1863     {
1864     Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");
1865     fclose(fileHandle_p);
1866 gross 1028 TMPMEMFREE(isCellCentered);
1867 dhawcroft 793 return;
1868     }
1869     nCompReqd = 3;
1870     }
1871     else
1872     {
1873     shape=getDataPointShape(data_pp[i_data], 0);
1874     if (shape>3 || shape != getDataPointShape(data_pp[i_data], 1))
1875     {
1876     Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
1877     fclose(fileHandle_p);
1878 gross 1028 TMPMEMFREE(isCellCentered);
1879 dhawcroft 793 return;
1880     }
1881     nCompReqd = 9;
1882     }
1883 gross 850 fprintf(fileHandle_p, "<DataArray Name=\"%s\" type=\"Float64\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);
1884 dhawcroft 793 /* write out the data */
1885     /* if the number of required components is more than the number
1886     * of actual components, pad with zeros
1887     */
1888 gross 1028 do_write=TRUE;
1889 dhawcroft 793 for (i=0; i<mesh_p->Nodes->numNodes; i++)
1890     {
1891     if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1892     {
1893     if (mesh_p->Nodes->toReduced[i]>=0)
1894     {
1895     switch(getFunctionSpaceType(data_pp[i_data]))
1896     {
1897     case FINLEY_DEGREES_OF_FREEDOM:
1898     values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);
1899     break;
1900     case FINLEY_REDUCED_DEGREES_OF_FREEDOM:
1901     values = getSampleData(data_pp[i_data],mesh_p->Nodes->reducedDegreeOfFreedom[i]);
1902     break;
1903     case FINLEY_NODES:
1904     values = getSampleData(data_pp[i_data],i);
1905     break;
1906     }
1907     do_write=TRUE;
1908     }
1909     else
1910     {
1911     do_write=FALSE;
1912     }
1913     }
1914     else
1915     {
1916     do_write=TRUE;
1917     switch(getFunctionSpaceType(data_pp[i_data]))
1918     {
1919     case FINLEY_DEGREES_OF_FREEDOM:
1920     values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);
1921     break;
1922     case FINLEY_NODES:
1923     values = getSampleData(data_pp[i_data],i);
1924     break;
1925     }
1926     }
1927     /* write the data different ways for scalar, vector and tensor */
1928     if (do_write)
1929     {
1930     if (nCompReqd == 1)
1931     {
1932     fprintf(fileHandle_p, " %e", values[0]);
1933     }
1934     else if (nCompReqd == 3)
1935     {
1936     for (m=0; m<shape; m++) fprintf(fileHandle_p, " %e", values[m]);
1937     for (m=0; m<nCompReqd-shape; m++) fprintf(fileHandle_p, " %e", 0.);
1938     }
1939     else if (nCompReqd == 9)
1940     {
1941     /* tensor data, so have a 3x3 matrix to output as a row
1942     * of 9 data points */
1943     count = 0;
1944     for (m=0; m<shape; m++)
1945     {
1946     for (n=0; n<shape; n++)
1947     {
1948     fprintf(fileHandle_p, " %e", values[count]);
1949     count++;
1950 jgs 153 }
1951 dhawcroft 793 for (n=0; n<3-shape; n++) fprintf(fileHandle_p, " %e", 0.);
1952     }
1953     for (m=0; m<3-shape; m++)
1954     for (n=0; n<3; n++) fprintf(fileHandle_p, " %e", 0.);
1955 jgs 153 }
1956 dhawcroft 793 fprintf(fileHandle_p, "\n");
1957     }
1958     }
1959     fprintf(fileHandle_p, "</DataArray>\n");
1960     }
1961     }
1962     fprintf(fileHandle_p, "</PointData>\n");
1963 jgs 153 }
1964 jgs 113 /* finish off the piece */
1965     fprintf(fileHandle_p, "</Piece>\n");
1966 jgs 110
1967     fprintf(fileHandle_p, "</UnstructuredGrid>\n");
1968     /* write the xml footer */
1969     fprintf(fileHandle_p, "</VTKFile>\n");
1970     /* close the file */
1971     fclose(fileHandle_p);
1972 gross 1028 TMPMEMFREE(isCellCentered);
1973 jgs 110 return;
1974     }
1975 dhawcroft 793 #endif
1976    

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26