/[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 850 - (hide annotations)
Sun Sep 17 23:27:00 2006 UTC (12 years, 10 months ago) by gross
File MIME type: text/plain
File size: 60607 byte(s)
some problems with vtk writer fixed
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 dhawcroft 818 rtmp = 0.;
870     for (j=0; j<numPointsPerSample; j++) rtmp += values[INDEX2(n,j,nComp)];
871     sampleAvg[k] = rtmp/numPointsPerSample;
872     }
873     // if the number of required components is more than the number
874     // of actual components, pad with zeros
875    
876     // probably only need to get shape of first element
877     // write the data different ways for scalar, vector and tensor
878     if (nCompReqd == 1)
879     {
880     sprintf(tmpbuf, " %e", sampleAvg[0]);
881     __STRCAT(largebuf,tmpbuf,tsz)
882     }
883     else if (nCompReqd == 3)
884     {
885     // write out the data
886     for (m=0; m<shape; m++)
887 dhawcroft 793 {
888 dhawcroft 818 sprintf(tmpbuf, " %e", sampleAvg[m]);
889     __STRCAT(largebuf,tmpbuf,tsz)
890 dhawcroft 793 }
891 dhawcroft 818 for (m=0; m<nCompReqd-shape; m++)
892 dhawcroft 793 {
893 dhawcroft 818 __STRCAT(largebuf,zero,tsz)
894 dhawcroft 793 }
895     }
896 dhawcroft 818 else if (nCompReqd == 9)
897 dhawcroft 793 {
898 dhawcroft 818 // tensor data, so have a 3x3 matrix to output as a row
899     // of 9 data points
900     count = 0;
901     for (m=0; m<shape; m++)
902 dhawcroft 793 {
903 dhawcroft 818 for (n=0; n<shape; n++)
904 dhawcroft 793 {
905 dhawcroft 818 sprintf(tmpbuf, " %e", sampleAvg[count]);
906     __STRCAT(largebuf,tmpbuf,tsz)
907     count++;
908 dhawcroft 793 }
909 dhawcroft 818 for (n=0; n<3-shape; n++)
910 dhawcroft 793 {
911 dhawcroft 818 __STRCAT(largebuf,zero,tsz)
912 dhawcroft 793 }
913     }
914 dhawcroft 818 for (m=0; m<3-shape; m++)
915     for (n=0; n<3; n++)
916 dhawcroft 793 {
917 dhawcroft 818 __STRCAT(largebuf,zero,tsz)
918 dhawcroft 793 }
919     }
920 dhawcroft 818 __STRCAT(largebuf,newline,tsz)
921 dhawcroft 793 }
922 dhawcroft 818 largebuf[tsz] = '\0';
923 dhawcroft 793 MPI_File_write_ordered(fh,largebuf,tsz,MPI_CHAR,&status);
924     MEMFREE(largebuf);
925     if( myRank == 0)
926     {
927     char *tag = "</DataArray>\n";
928     MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);
929     MPI_Wait(&req,&status);
930     }
931 dhawcroft 818
932 dhawcroft 793 }
933     }
934 dhawcroft 818 // closing celldata tag
935 dhawcroft 793 if(myRank == 0)
936     {
937 dhawcroft 818 char* tag = "</CellData>\n";
938 dhawcroft 793 MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);
939     MPI_Wait(&req,&status);
940     }
941 dhawcroft 818
942     MPIO_DEBUG(" Done Writing Cell Data ")
943 dhawcroft 793 }
944    
945 dhawcroft 818
946     // Write Point Data Header Tags
947     if( myRank == 0)
948 dhawcroft 793 {
949 dhawcroft 818 char header[600];
950     char tmpBuf[50];
951    
952     if (write_pointdata)
953 dhawcroft 793 {
954 dhawcroft 818 MPIO_DEBUG(" Writing Pointdata... ")
955 dhawcroft 793 // mark the active data arrays
956     bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
957 dhawcroft 818 sprintf(header, "<PointData");
958 dhawcroft 793 for (i_data =0 ;i_data<num_data;++i_data)
959     {
960 dhawcroft 818 if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])
961 dhawcroft 793 {
962     // if the rank == 0: --> scalar data
963     // if the rank == 1: --> vector data
964     // if the rank == 2: --> tensor data
965    
966     switch(getDataPointRank(data_pp[i_data]))
967     {
968     case 0:
969     if (! set_scalar)
970     {
971     sprintf(tmpBuf," Scalars=\"%s\"",names_p[i_data]);
972     strcat(header,tmpBuf);
973     set_scalar=TRUE;
974     }
975     break;
976     case 1:
977     if (! set_vector)
978     {
979     sprintf(tmpBuf," Vectors=\"%s\"",names_p[i_data]);
980 dhawcroft 818 strcat(header,tmpBuf);
981 dhawcroft 793 set_vector=TRUE;
982     }
983     break;
984     case 2:
985     if (! set_tensor)
986     {
987     sprintf(tmpBuf," Tensors=\"%s\"",names_p[i_data]);
988 dhawcroft 818 strcat(header,tmpBuf);
989 dhawcroft 793 set_tensor=TRUE;
990     }
991     break;
992     default:
993     sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);
994     Finley_setError(VALUE_ERROR,error_msg);
995     return;
996     }
997     }
998     }
999     strcat(header, ">\n");
1000     MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);
1001     MPI_Wait(&req,&status);
1002     }
1003     }
1004    
1005 dhawcroft 818 // write actual data
1006     if(write_pointdata)
1007 dhawcroft 793 {
1008     for (i_data =0 ;i_data<num_data;++i_data)
1009     {
1010 dhawcroft 818 if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])
1011 dhawcroft 793 {
1012     numPointsPerSample = elements->ReferenceElement->numQuadNodes;
1013     rank = getDataPointRank(data_pp[i_data]);
1014     nComp = getDataPointSize(data_pp[i_data]);
1015     nCompReqd=1; // the number of components required by vtk
1016     shape=0;
1017     if (rank == 0)
1018     {
1019     nCompReqd = 1;
1020     }
1021     else if (rank == 1)
1022     {
1023     shape=getDataPointShape(data_pp[i_data], 0);
1024     if (shape>3)
1025     {
1026     Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");
1027     return;
1028     }
1029     nCompReqd = 3;
1030     }
1031     else
1032     {
1033     shape=getDataPointShape(data_pp[i_data], 0);
1034     if (shape>3 || shape != getDataPointShape(data_pp[i_data], 1))
1035     {
1036     Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
1037     return;
1038     }
1039     nCompReqd = 9;
1040     }
1041    
1042     if( myRank == 0)
1043     {
1044     char header[250];
1045 gross 850 sprintf(header, "<DataArray Name=\"%s\" type=\"Float64\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);
1046 dhawcroft 793 MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);
1047     MPI_Wait(&req,&status);
1048     }
1049 dhawcroft 818 // write out the data
1050     // if the number of required components is more than the number
1051     // of actual components, pad with zeros
1052 dhawcroft 793
1053 dhawcroft 818 char tmpbuf[15];
1054     char* largebuf = MEMALLOC(nCompReqd*numLocalNodes*15 + numLocal*nCompReqd + nCompReqd + 15,char);
1055 dhawcroft 793 largebuf[0] = '\0';
1056 dhawcroft 818 bool_t do_write=TRUE;
1057 dhawcroft 793 size_t tsz = 0;
1058    
1059 dhawcroft 818 for(k=0;k < nodeCache.size;k++)
1060 dhawcroft 793 {
1061 dhawcroft 818 i = nodeCache.values[k];
1062 dhawcroft 793
1063 dhawcroft 818 if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1064 dhawcroft 793 {
1065 dhawcroft 818 if (mesh_p->Nodes->toReduced[i]>=0)
1066     {
1067     switch(getFunctionSpaceType(data_pp[i_data]))
1068     {
1069     case FINLEY_DEGREES_OF_FREEDOM:
1070     values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);
1071     break;
1072     case FINLEY_REDUCED_DEGREES_OF_FREEDOM:
1073     values = getSampleData(data_pp[i_data],mesh_p->Nodes->reducedDegreeOfFreedom[i]);
1074     break;
1075     case FINLEY_NODES:
1076     values = getSampleData(data_pp[i_data],i);
1077     break;
1078     }
1079     do_write=TRUE;
1080     }
1081     else
1082     {
1083     do_write=FALSE;
1084     }
1085 dhawcroft 793 }
1086 dhawcroft 818 else
1087 dhawcroft 793 {
1088 dhawcroft 818 do_write=TRUE;
1089     switch(getFunctionSpaceType(data_pp[i_data]))
1090     {
1091     case FINLEY_DEGREES_OF_FREEDOM:
1092     values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);
1093     break;
1094     case FINLEY_NODES:
1095     values = getSampleData(data_pp[i_data],i);
1096     break;
1097     }
1098 dhawcroft 793 }
1099 dhawcroft 818 // write the data different ways for scalar, vector and tensor
1100     if (do_write)
1101 dhawcroft 793 {
1102 dhawcroft 818 if (nCompReqd == 1)
1103 dhawcroft 793 {
1104 dhawcroft 818 sprintf(tmpbuf," %e", values[0]);
1105     __STRCAT(largebuf,tmpbuf,tsz)
1106 dhawcroft 793 }
1107 dhawcroft 818 else if (nCompReqd == 3)
1108 dhawcroft 793 {
1109 dhawcroft 818 for (m=0; m<shape; m++)
1110     {
1111 dhawcroft 793
1112 dhawcroft 818 sprintf(tmpbuf," %e",values[m]);
1113     __STRCAT(largebuf,tmpbuf,tsz)
1114     }
1115     for (m=0; m<nCompReqd-shape; m++)
1116     {
1117     __STRCAT(largebuf,zero,tsz)
1118     }
1119 dhawcroft 793 }
1120 dhawcroft 818 else if (nCompReqd == 9)
1121 dhawcroft 793 {
1122 dhawcroft 818 // tensor data, so have a 3x3 matrix to output as a row
1123     // of 9 data points
1124     count = 0;
1125     for (m=0; m<shape; m++)
1126 dhawcroft 793 {
1127 dhawcroft 818 for (n=0; n<shape; n++)
1128     {
1129     sprintf(tmpbuf," %e", values[count]);
1130     __STRCAT(largebuf,tmpbuf,tsz)
1131     count++;
1132     }
1133     for (n=0; n<3-shape; n++)
1134     {
1135     __STRCAT(largebuf,zero,tsz)
1136     }
1137 dhawcroft 793 }
1138 dhawcroft 818 for (m=0; m<3-shape; m++)
1139 dhawcroft 793 {
1140 dhawcroft 818 for (n=0; n<3; n++)
1141     {
1142     __STRCAT(largebuf,zero,tsz)
1143     }
1144 dhawcroft 793 }
1145     }
1146 dhawcroft 818 __STRCAT(largebuf,newline,tsz)
1147 dhawcroft 793 }
1148 dhawcroft 818
1149 dhawcroft 793 }
1150 dhawcroft 818 // Write out local data
1151    
1152     largebuf[tsz] = '\0';
1153 dhawcroft 793 MPI_File_write_ordered(fh,largebuf,tsz,MPI_CHAR,&status);
1154     MEMFREE(largebuf);
1155     if( myRank == 0)
1156     {
1157     char *tag = "</DataArray>\n";
1158     MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);
1159     MPI_Wait(&req,&status);
1160     }
1161     }
1162     }
1163 dhawcroft 818 // Finish off with closing tag
1164 dhawcroft 793 if(myRank == 0)
1165     {
1166 dhawcroft 818 char* tag = "</PointData>\n";
1167 dhawcroft 793 MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);
1168     MPI_Wait(&req,&status);
1169     }
1170 dhawcroft 818 MPIO_DEBUG(" Done Writing Pointdata ")
1171 dhawcroft 793 }
1172 dhawcroft 818 // end write_pointdata
1173 dhawcroft 793
1174 dhawcroft 818 // tag and bag...
1175 dhawcroft 793 if (myRank == 0)
1176     {
1177     char *footer = "</Piece>\n</UnstructuredGrid>\n</VTKFile>";
1178     MPI_File_iwrite_shared(fh,footer,strlen(footer),MPI_CHAR,&req);
1179     MPI_Wait(&req,&status);
1180     }
1181    
1182     MEMFREE(nodesGlobal);
1183     MEMFREE(nodeCache.values);
1184     MEMFREE(elementCache.values);
1185     #ifdef MPIO_HINTS
1186     MPI_Info_free(&infoHints);
1187 dhawcroft 796 #undef MPIO_HINTS
1188 dhawcroft 793 #endif
1189     MPI_File_close(&fh);
1190     MPIO_DEBUG(" ***** Exit saveVTK ***** ")
1191 dhawcroft 796 #undef __STRCAT
1192 dhawcroft 793 }
1193    
1194     #undef MPIO_DEBUG
1195     #else
1196    
1197    
1198    
1199    
1200     void Finley_Mesh_saveVTK(const char * filename_p, Finley_Mesh *mesh_p, const dim_t num_data,char* *names_p, escriptDataC* *data_pp)
1201     {
1202     char error_msg[LenErrorMsg_MAX];
1203 jgs 110 /* if there is no mesh we just return */
1204     if (mesh_p==NULL) return;
1205 jgs 153
1206 gross 687 int i, j, k, numVTKNodesPerElement,i_data,m, count, n, rank,shape, numPoints, cellType, numCells,
1207 dhawcroft 793 nDim, numPointsPerSample, nComp, nCompReqd;
1208 gross 687
1209 jgs 150 index_t j2;
1210 jgs 113 double* values, rtmp;
1211 gross 687 char elemTypeStr[32];
1212 dhawcroft 793
1213 jgs 153 /* open the file and check handle */
1214    
1215     FILE * fileHandle_p = fopen(filename_p, "w");
1216 dhawcroft 793 if (fileHandle_p==NULL)
1217     {
1218 jgs 153 sprintf(error_msg, "saveVTK: File %s could not be opened for writing.", filename_p);
1219     Finley_setError(IO_ERROR,error_msg);
1220     return;
1221 jgs 150 }
1222 jgs 153 /* find the mesh type to be written */
1223     int nodetype=FINLEY_DEGREES_OF_FREEDOM;
1224     int elementtype=FINLEY_UNKNOWN;
1225     bool_t isCellCentered[num_data],write_celldata=FALSE,write_pointdata=FALSE;
1226 dhawcroft 793 for (i_data=0;i_data<num_data;++i_data)
1227     {
1228     if (! isEmpty(data_pp[i_data]))
1229     {
1230     switch(getFunctionSpaceType(data_pp[i_data]))
1231     {
1232     case FINLEY_DEGREES_OF_FREEDOM:
1233     nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1234     if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)
1235     {
1236     elementtype=FINLEY_ELEMENTS;
1237 jgs 153 }
1238 dhawcroft 793 else
1239     {
1240     Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1241     fclose(fileHandle_p);
1242     return;
1243 jgs 153 }
1244 dhawcroft 793 isCellCentered[i_data]=FALSE;
1245     break;
1246     case FINLEY_REDUCED_DEGREES_OF_FREEDOM:
1247     nodetype = FINLEY_REDUCED_DEGREES_OF_FREEDOM;
1248     if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)
1249     {
1250     elementtype=FINLEY_ELEMENTS;
1251     }
1252     else
1253     {
1254     Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1255     fclose(fileHandle_p);
1256     return;
1257     }
1258     isCellCentered[i_data]=FALSE;
1259     break;
1260     case FINLEY_NODES:
1261     nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1262     if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)
1263     {
1264     elementtype=FINLEY_ELEMENTS;
1265     }
1266     else
1267     {
1268     Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1269     fclose(fileHandle_p);
1270     return;
1271     }
1272     isCellCentered[i_data]=FALSE;
1273     break;
1274     case FINLEY_ELEMENTS:
1275     nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1276     if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)
1277     {
1278     elementtype=FINLEY_ELEMENTS;
1279     }
1280     else
1281     {
1282     Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1283     fclose(fileHandle_p);
1284     return;
1285     }
1286     isCellCentered[i_data]=TRUE;
1287     break;
1288     case FINLEY_FACE_ELEMENTS:
1289     nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1290     if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_FACE_ELEMENTS)
1291     {
1292     elementtype=FINLEY_FACE_ELEMENTS;
1293     }
1294     else
1295     {
1296     Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1297     fclose(fileHandle_p);
1298     return;
1299     }
1300     isCellCentered[i_data]=TRUE;
1301     break;
1302     case FINLEY_POINTS:
1303     nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1304     if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_POINTS)
1305     {
1306     elementtype=FINLEY_POINTS;
1307     }
1308     else
1309     {
1310     Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1311     fclose(fileHandle_p);
1312     return;
1313     }
1314     isCellCentered[i_data]=TRUE;
1315     break;
1316     case FINLEY_CONTACT_ELEMENTS_1:
1317     nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1318     if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1)
1319     {
1320     elementtype=FINLEY_CONTACT_ELEMENTS_1;
1321     }
1322     else
1323     {
1324     Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1325     fclose(fileHandle_p);
1326     return;
1327     }
1328     isCellCentered[i_data]=TRUE;
1329     break;
1330     case FINLEY_CONTACT_ELEMENTS_2:
1331     nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1332     if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1)
1333     {
1334     elementtype=FINLEY_CONTACT_ELEMENTS_1;
1335     }
1336     else
1337     {
1338     Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1339     fclose(fileHandle_p);
1340     return;
1341     }
1342     isCellCentered[i_data]=TRUE;
1343     break;
1344     default:
1345     sprintf(error_msg,"saveVTK: Finley does not know anything about function space type %d",getFunctionSpaceType(data_pp[i_data]));
1346     Finley_setError(TYPE_ERROR,error_msg);
1347     fclose(fileHandle_p);
1348     return;
1349     }
1350     if (isCellCentered[i_data])
1351     {
1352     write_celldata=TRUE;
1353     }
1354     else
1355     {
1356     write_pointdata=TRUE;
1357     }
1358     }
1359 jgs 153 }
1360     /* select nomber of points and the mesh component */
1361 gross 687 numPoints = mesh_p->Nodes->numNodes;
1362 dhawcroft 793 if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1363     {
1364     numPoints = mesh_p->Nodes->reducedNumNodes;
1365 jgs 153 }
1366 dhawcroft 793 else
1367     {
1368     numPoints = mesh_p->Nodes->numNodes;
1369     }
1370 jgs 153 if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;
1371     Finley_ElementFile* elements=NULL;
1372 dhawcroft 793 switch(elementtype)
1373     {
1374     case FINLEY_ELEMENTS:
1375     elements=mesh_p->Elements;
1376     break;
1377     case FINLEY_FACE_ELEMENTS:
1378     elements=mesh_p->FaceElements;
1379     break;
1380     case FINLEY_POINTS:
1381     elements=mesh_p->Points;
1382     break;
1383     case FINLEY_CONTACT_ELEMENTS_1:
1384     elements=mesh_p->ContactElements;
1385     break;
1386 jgs 113 }
1387 dhawcroft 793 if (elements==NULL)
1388     {
1389     Finley_setError(SYSTEM_ERROR,"saveVTK: undefined element file");
1390     fclose(fileHandle_p);
1391     return;
1392 jgs 113 }
1393 jgs 153 /* map finley element type to VTK element type */
1394 dhawcroft 793 numCells = elements->numElements;
1395 jgs 153 ElementTypeId TypeId;
1396 dhawcroft 793 if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1397     {
1398     TypeId = elements->LinearReferenceElement->Type->TypeId;
1399 jgs 113 }
1400 dhawcroft 793 else
1401     {
1402     TypeId = elements->ReferenceElement->Type->TypeId;
1403     }
1404 jgs 110
1405 dhawcroft 793 switch(TypeId)
1406     {
1407     case Point1:
1408     case Line2Face:
1409     case Line3Face:
1410     case Point1_Contact:
1411     case Line2Face_Contact:
1412     case Line3Face_Contact:
1413     cellType = VTK_VERTEX;
1414     numVTKNodesPerElement = 1;
1415     strcpy(elemTypeStr, "VTK_VERTEX");
1416     break;
1417 jgs 153
1418 dhawcroft 793 case Line2:
1419     case Tri3Face:
1420     case Rec4Face:
1421     case Line2_Contact:
1422     case Tri3_Contact:
1423     case Tri3Face_Contact:
1424     case Rec4Face_Contact:
1425     cellType = VTK_LINE;
1426     numVTKNodesPerElement = 2;
1427     strcpy(elemTypeStr, "VTK_LINE");
1428     break;
1429 jgs 153
1430 dhawcroft 793 case Tri3:
1431     case Tet4Face:
1432     case Tet4Face_Contact:
1433     cellType = VTK_TRIANGLE;
1434     numVTKNodesPerElement = 3;
1435     strcpy(elemTypeStr, "VTK_TRIANGLE");
1436     break;
1437 jgs 153
1438 dhawcroft 793 case Rec4:
1439     case Hex8Face:
1440     case Rec4_Contact:
1441     case Hex8Face_Contact:
1442     cellType = VTK_QUAD;
1443     numVTKNodesPerElement = 4;
1444     strcpy(elemTypeStr, "VTK_QUAD");
1445     break;
1446 jgs 153
1447 dhawcroft 793 case Tet4:
1448     cellType = VTK_TETRA;
1449     numVTKNodesPerElement = 4;
1450     strcpy(elemTypeStr, "VTK_TETRA");
1451     break;
1452 jgs 153
1453 dhawcroft 793 case Hex8:
1454     cellType = VTK_HEXAHEDRON;
1455     numVTKNodesPerElement = 8;
1456     strcpy(elemTypeStr, "VTK_HEXAHEDRON");
1457     break;
1458 jgs 153
1459 dhawcroft 793 case Line3:
1460     case Tri6Face:
1461     case Rec8Face:
1462     case Line3_Contact:
1463     case Tri6Face_Contact:
1464     case Rec8Face_Contact:
1465     cellType = VTK_QUADRATIC_EDGE;
1466     numVTKNodesPerElement = 3;
1467     strcpy(elemTypeStr, "VTK_QUADRATIC_EDGE");
1468     break;
1469 jgs 153
1470 dhawcroft 793 case Tri6:
1471     case Tet10Face:
1472     case Tri6_Contact:
1473     case Tet10Face_Contact:
1474     cellType = VTK_QUADRATIC_TRIANGLE;
1475     numVTKNodesPerElement = 6;
1476     strcpy(elemTypeStr, "VTK_QUADRATIC_TRIANGLE");
1477     break;
1478 jgs 153
1479 dhawcroft 793 case Rec8:
1480     case Hex20Face:
1481     case Rec8_Contact:
1482     case Hex20Face_Contact:
1483     cellType = VTK_QUADRATIC_QUAD;
1484     numVTKNodesPerElement = 8;
1485     strcpy(elemTypeStr, "VTK_QUADRATIC_QUAD");
1486     break;
1487 jgs 153
1488 dhawcroft 793 case Tet10:
1489     cellType = VTK_QUADRATIC_TETRA;
1490     numVTKNodesPerElement = 10;
1491     strcpy(elemTypeStr, "VTK_QUADRATIC_TETRA");
1492     break;
1493 jgs 153
1494 dhawcroft 793 case Hex20:
1495     cellType = VTK_QUADRATIC_HEXAHEDRON;
1496     numVTKNodesPerElement = 20;
1497     strcpy(elemTypeStr, "VTK_QUADRATIC_HEXAHEDRON");
1498     break;
1499 jgs 153
1500 dhawcroft 793 default:
1501     sprintf(error_msg, "saveVTK: Element type %s is not supported by VTK",elements->ReferenceElement->Type->Name);
1502     Finley_setError(VALUE_ERROR,error_msg);
1503     fclose(fileHandle_p);
1504     return;
1505     }
1506 jgs 153 /* xml header */
1507     fprintf(fileHandle_p, "<?xml version=\"1.0\"?>\n");
1508     fprintf(fileHandle_p, "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\">\n");
1509 jgs 110
1510 jgs 153 /* finley uses an unstructured mesh, so UnstructuredGrid *should* work */
1511     fprintf(fileHandle_p, "<UnstructuredGrid>\n");
1512    
1513     /* is there only one "piece" to the data?? */
1514     fprintf(fileHandle_p, "<Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n",numPoints, numCells);
1515     /* now for the points; equivalent to positions section in saveDX() */
1516     /* "The points element explicitly defines coordinates for each point
1517 dhawcroft 793 * individually. It contains one DataArray element describing an array
1518     * with three components per value, each specifying the coordinates of one
1519     * point" - from Vtk User's Guide
1520     */
1521 jgs 153 fprintf(fileHandle_p, "<Points>\n");
1522 dhawcroft 793 /*
1523     * the reason for this if statement is explained in the long comment below
1524     */
1525 gross 687 nDim = mesh_p->Nodes->numDim;
1526 gross 850 fprintf(fileHandle_p, "<DataArray NumberOfComponents=\"%d\" type=\"Float64\" format=\"ascii\">\n",MAX(3,nDim));
1527 jgs 153 /* vtk/mayavi doesn't like 2D data, it likes 3D data with a degenerate
1528 dhawcroft 793 * third dimension to handle 2D data (like a sheet of paper). So, if
1529     * nDim is 2, we have to append zeros to the array to get this third
1530     * dimension, and keep the visualisers happy.
1531     * Indeed, if nDim is less than 3, must pad all empty dimensions, so
1532     * that the total number of dims is 3.
1533 jgs 153 */
1534 dhawcroft 793 if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1535     {
1536     for (i = 0; i < mesh_p->Nodes->numNodes; i++)
1537     {
1538     if (mesh_p->Nodes->toReduced[i]>=0)
1539     {
1540     for (j = 0; j < nDim; j++) fprintf(fileHandle_p, " %e",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);
1541     for (k=0; k<3-nDim; k++) fprintf(fileHandle_p, " %e",0.);
1542     fprintf(fileHandle_p, "\n");
1543     }
1544     }
1545     }
1546     else
1547     {
1548     for (i = 0; i < mesh_p->Nodes->numNodes; i++)
1549     {
1550    
1551     for (j = 0; j < nDim; j++) fprintf(fileHandle_p, " %e",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);
1552     for (k=0; k<3-nDim; k++) fprintf(fileHandle_p, " %e",0.);
1553     fprintf(fileHandle_p, "\n");
1554     }
1555    
1556     }
1557 jgs 153 fprintf(fileHandle_p, "</DataArray>\n");
1558     fprintf(fileHandle_p, "</Points>\n");
1559 jgs 113
1560 jgs 153 /* write out the DataArray element for the connectivity */
1561 jgs 113
1562 jgs 153 int NN = elements->ReferenceElement->Type->numNodes;
1563     fprintf(fileHandle_p, "<Cells>\n");
1564     fprintf(fileHandle_p, "<DataArray Name=\"connectivity\" type=\"Int32\" format=\"ascii\">\n");
1565 jgs 113
1566 dhawcroft 793 if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1567     {
1568     for (i = 0; i < numCells; i++)
1569     {
1570     for (j = 0; j < numVTKNodesPerElement; j++)
1571     fprintf(fileHandle_p,"%d ",mesh_p->Nodes->toReduced[elements->Nodes[INDEX2(elements->ReferenceElement->Type->linearNodes[j], i, NN)]]);
1572     fprintf(fileHandle_p, "\n");
1573     }
1574     }
1575     else if (VTK_QUADRATIC_HEXAHEDRON==cellType)
1576     {
1577     for (i = 0; i < numCells; i++)
1578     {
1579     fprintf(fileHandle_p,"%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d\n",
1580     elements->Nodes[INDEX2(0, i, NN)],
1581     elements->Nodes[INDEX2(1, i, NN)],
1582     elements->Nodes[INDEX2(2, i, NN)],
1583     elements->Nodes[INDEX2(3, i, NN)],
1584     elements->Nodes[INDEX2(4, i, NN)],
1585     elements->Nodes[INDEX2(5, i, NN)],
1586     elements->Nodes[INDEX2(6, i, NN)],
1587     elements->Nodes[INDEX2(7, i, NN)],
1588     elements->Nodes[INDEX2(8, i, NN)],
1589     elements->Nodes[INDEX2(9, i, NN)],
1590     elements->Nodes[INDEX2(10, i, NN)],
1591     elements->Nodes[INDEX2(11, i, NN)],
1592     elements->Nodes[INDEX2(16, i, NN)],
1593     elements->Nodes[INDEX2(17, i, NN)],
1594     elements->Nodes[INDEX2(18, i, NN)],
1595     elements->Nodes[INDEX2(19, i, NN)],
1596     elements->Nodes[INDEX2(12, i, NN)],
1597     elements->Nodes[INDEX2(13, i, NN)],
1598     elements->Nodes[INDEX2(14, i, NN)],
1599     elements->Nodes[INDEX2(15, i, NN)]);
1600     }
1601     }
1602     else if (numVTKNodesPerElement!=NN)
1603     {
1604     for (i = 0; i < numCells; i++)
1605     {
1606     for (j = 0; j < numVTKNodesPerElement; j++) fprintf(fileHandle_p,"%d ", elements->Nodes[INDEX2(elements->ReferenceElement->Type->geoNodes[j], i, NN)]);
1607     fprintf(fileHandle_p, "\n");
1608     }
1609     }
1610     else
1611     {
1612     for (i = 0; i < numCells; i++)
1613     {
1614     for (j = 0; j < numVTKNodesPerElement; j++) fprintf(fileHandle_p,"%d ", elements->Nodes[INDEX2(j, i, NN)]);
1615     fprintf(fileHandle_p, "\n");
1616     }
1617     }
1618 jgs 153 fprintf(fileHandle_p, "</DataArray>\n");
1619 jgs 110
1620 jgs 153 /* write out the DataArray element for the offsets */
1621     fprintf(fileHandle_p, "<DataArray Name=\"offsets\" type=\"Int32\" format=\"ascii\">\n");
1622     for (i=numVTKNodesPerElement; i<=numCells*numVTKNodesPerElement; i+=numVTKNodesPerElement) fprintf(fileHandle_p, "%d\n", i);
1623     fprintf(fileHandle_p, "</DataArray>\n");
1624    
1625     /* write out the DataArray element for the types */
1626     fprintf(fileHandle_p, "<DataArray Name=\"types\" type=\"UInt8\" format=\"ascii\">\n");
1627     for (i=0; i<numCells; i++) fprintf(fileHandle_p, "%d\n", cellType);
1628     fprintf(fileHandle_p, "</DataArray>\n");
1629    
1630     /* finish off the <Cells> element */
1631     fprintf(fileHandle_p, "</Cells>\n");
1632    
1633     /* cell data */
1634 dhawcroft 793 if (write_celldata)
1635     {
1636     /* mark the active data arrays */
1637     bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
1638     fprintf(fileHandle_p, "<CellData");
1639     for (i_data =0 ;i_data<num_data;++i_data)
1640     {
1641     if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data])
1642     {
1643     /* if the rank == 0: --> scalar data
1644     * if the rank == 1: --> vector data
1645     * if the rank == 2: --> tensor data
1646     */
1647     switch(getDataPointRank(data_pp[i_data]))
1648     {
1649     case 0:
1650     if (! set_scalar)
1651     {
1652     fprintf(fileHandle_p," Scalars=\"%s\"",names_p[i_data]);
1653     set_scalar=TRUE;
1654     }
1655     break;
1656     case 1:
1657     if (! set_vector)
1658     {
1659     fprintf(fileHandle_p," Vectors=\"%s\"",names_p[i_data]);
1660     set_vector=TRUE;
1661     }
1662     break;
1663     case 2:
1664     if (! set_tensor)
1665     {
1666     fprintf(fileHandle_p," Tensors=\"%s\"",names_p[i_data]);
1667     set_tensor=TRUE;
1668     }
1669     break;
1670     default:
1671     sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);
1672     Finley_setError(VALUE_ERROR,error_msg);
1673     fclose(fileHandle_p);
1674     return;
1675     }
1676     }
1677     }
1678     fprintf(fileHandle_p, ">\n");
1679     /* write the arrays */
1680     for (i_data =0 ;i_data<num_data;++i_data)
1681     {
1682     if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data])
1683     {
1684     numPointsPerSample = elements->ReferenceElement->numQuadNodes;
1685     rank = getDataPointRank(data_pp[i_data]);
1686     nComp = getDataPointSize(data_pp[i_data]);
1687     nCompReqd=1; /* the number of components required by vtk */
1688     shape=0;
1689     if (rank == 0)
1690     {
1691     nCompReqd = 1;
1692     }
1693     else if (rank == 1)
1694     {
1695     shape=getDataPointShape(data_pp[i_data], 0);
1696     if (shape>3)
1697     {
1698     Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");
1699     fclose(fileHandle_p);
1700     return;
1701     }
1702     nCompReqd = 3;
1703     }
1704     else
1705     {
1706     shape=getDataPointShape(data_pp[i_data], 0);
1707     if (shape>3 || shape != getDataPointShape(data_pp[i_data], 1))
1708     {
1709     Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
1710     fclose(fileHandle_p);
1711     return;
1712     }
1713     nCompReqd = 9;
1714     }
1715 gross 850 fprintf(fileHandle_p, "<DataArray Name=\"%s\" type=\"Float64\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);
1716 dhawcroft 793
1717     double sampleAvg[nComp];
1718     for (i=0; i<numCells; i++)
1719     {
1720     values = getSampleData(data_pp[i_data], i);
1721     /* averaging over the number of points in the sample */
1722     for (k=0; k<nComp; k++)
1723     {
1724     rtmp = 0.;
1725     for (j=0; j<numPointsPerSample; j++) rtmp += values[INDEX2(k,j,nComp)];
1726     sampleAvg[k] = rtmp/numPointsPerSample;
1727     }
1728     /* if the number of required components is more than the number
1729     * of actual components, pad with zeros
1730     */
1731     /* probably only need to get shape of first element */
1732     /* write the data different ways for scalar, vector and tensor */
1733     if (nCompReqd == 1)
1734     {
1735     fprintf(fileHandle_p, " %e", sampleAvg[0]);
1736     }
1737     else if (nCompReqd == 3)
1738     {
1739     /* write out the data */
1740     for (m=0; m<shape; m++) fprintf(fileHandle_p, " %e", sampleAvg[m]);
1741     for (m=0; m<nCompReqd-shape; m++) fprintf(fileHandle_p, " %e", 0.);
1742     }
1743     else if (nCompReqd == 9)
1744     {
1745     /* tensor data, so have a 3x3 matrix to output as a row
1746     * of 9 data points */
1747     count = 0;
1748     for (m=0; m<shape; m++)
1749     {
1750     for (n=0; n<shape; n++)
1751     {
1752     fprintf(fileHandle_p, " %e", sampleAvg[count]);
1753     count++;
1754     }
1755     for (n=0; n<3-shape; n++) fprintf(fileHandle_p, " %e", 0.);
1756 jgs 150 }
1757 dhawcroft 793 for (m=0; m<3-shape; m++)
1758     for (n=0; n<3; n++) fprintf(fileHandle_p, " %e", 0.);
1759     }
1760     fprintf(fileHandle_p, "\n");
1761     }
1762     fprintf(fileHandle_p, "</DataArray>\n");
1763     }
1764     }
1765     fprintf(fileHandle_p, "</CellData>\n");
1766 jgs 110 }
1767 jgs 153 /* point data */
1768 dhawcroft 793 if (write_pointdata)
1769     {
1770     /* mark the active data arrays */
1771     bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
1772     fprintf(fileHandle_p, "<PointData");
1773     for (i_data =0 ;i_data<num_data;++i_data)
1774     {
1775     if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])
1776     {
1777     /* if the rank == 0: --> scalar data
1778     * if the rank == 1: --> vector data
1779     * if the rank == 2: --> tensor data
1780     */
1781     switch(getDataPointRank(data_pp[i_data]))
1782     {
1783     case 0:
1784     if (! set_scalar)
1785     {
1786     fprintf(fileHandle_p," Scalars=\"%s\"",names_p[i_data]);
1787     set_scalar=TRUE;
1788     }
1789     break;
1790     case 1:
1791     if (! set_vector)
1792     {
1793     fprintf(fileHandle_p," Vectors=\"%s\"",names_p[i_data]);
1794     set_vector=TRUE;
1795     }
1796     break;
1797     case 2:
1798     if (! set_tensor)
1799     {
1800     fprintf(fileHandle_p," Tensors=\"%s\"",names_p[i_data]);
1801     set_tensor=TRUE;
1802     }
1803     break;
1804     default:
1805     sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);
1806     Finley_setError(VALUE_ERROR,error_msg);
1807     fclose(fileHandle_p);
1808     return;
1809     }
1810     }
1811     }
1812     fprintf(fileHandle_p, ">\n");
1813     /* write the arrays */
1814     for (i_data =0 ;i_data<num_data;++i_data)
1815     {
1816     if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])
1817     {
1818     numPointsPerSample = elements->ReferenceElement->numQuadNodes;
1819     rank = getDataPointRank(data_pp[i_data]);
1820     nComp = getDataPointSize(data_pp[i_data]);
1821     nCompReqd=1; /* the number of components required by vtk */
1822     shape=0;
1823     if (rank == 0)
1824     {
1825     nCompReqd = 1;
1826     }
1827     else if (rank == 1)
1828     {
1829     shape=getDataPointShape(data_pp[i_data], 0);
1830     if (shape>3)
1831     {
1832     Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");
1833     fclose(fileHandle_p);
1834     return;
1835     }
1836     nCompReqd = 3;
1837     }
1838     else
1839     {
1840     shape=getDataPointShape(data_pp[i_data], 0);
1841     if (shape>3 || shape != getDataPointShape(data_pp[i_data], 1))
1842     {
1843     Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
1844     fclose(fileHandle_p);
1845     return;
1846     }
1847     nCompReqd = 9;
1848     }
1849 gross 850 fprintf(fileHandle_p, "<DataArray Name=\"%s\" type=\"Float64\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);
1850 dhawcroft 793 /* write out the data */
1851     /* if the number of required components is more than the number
1852     * of actual components, pad with zeros
1853     */
1854     bool_t do_write=TRUE;
1855     for (i=0; i<mesh_p->Nodes->numNodes; i++)
1856     {
1857     if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1858     {
1859     if (mesh_p->Nodes->toReduced[i]>=0)
1860     {
1861     switch(getFunctionSpaceType(data_pp[i_data]))
1862     {
1863     case FINLEY_DEGREES_OF_FREEDOM:
1864     values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);
1865     break;
1866     case FINLEY_REDUCED_DEGREES_OF_FREEDOM:
1867     values = getSampleData(data_pp[i_data],mesh_p->Nodes->reducedDegreeOfFreedom[i]);
1868     break;
1869     case FINLEY_NODES:
1870     values = getSampleData(data_pp[i_data],i);
1871     break;
1872     }
1873     do_write=TRUE;
1874     }
1875     else
1876     {
1877     do_write=FALSE;
1878     }
1879     }
1880     else
1881     {
1882     do_write=TRUE;
1883     switch(getFunctionSpaceType(data_pp[i_data]))
1884     {
1885     case FINLEY_DEGREES_OF_FREEDOM:
1886     values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);
1887     break;
1888     case FINLEY_NODES:
1889     values = getSampleData(data_pp[i_data],i);
1890     break;
1891     }
1892     }
1893     /* write the data different ways for scalar, vector and tensor */
1894     if (do_write)
1895     {
1896     if (nCompReqd == 1)
1897     {
1898     fprintf(fileHandle_p, " %e", values[0]);
1899     }
1900     else if (nCompReqd == 3)
1901     {
1902     for (m=0; m<shape; m++) fprintf(fileHandle_p, " %e", values[m]);
1903     for (m=0; m<nCompReqd-shape; m++) fprintf(fileHandle_p, " %e", 0.);
1904     }
1905     else if (nCompReqd == 9)
1906     {
1907     /* tensor data, so have a 3x3 matrix to output as a row
1908     * of 9 data points */
1909     count = 0;
1910     for (m=0; m<shape; m++)
1911     {
1912     for (n=0; n<shape; n++)
1913     {
1914     fprintf(fileHandle_p, " %e", values[count]);
1915     count++;
1916 jgs 153 }
1917 dhawcroft 793 for (n=0; n<3-shape; n++) fprintf(fileHandle_p, " %e", 0.);
1918     }
1919     for (m=0; m<3-shape; m++)
1920     for (n=0; n<3; n++) fprintf(fileHandle_p, " %e", 0.);
1921 jgs 153 }
1922 dhawcroft 793 fprintf(fileHandle_p, "\n");
1923     }
1924     }
1925     fprintf(fileHandle_p, "</DataArray>\n");
1926     }
1927     }
1928     fprintf(fileHandle_p, "</PointData>\n");
1929 jgs 153 }
1930 jgs 113 /* finish off the piece */
1931     fprintf(fileHandle_p, "</Piece>\n");
1932 jgs 110
1933     fprintf(fileHandle_p, "</UnstructuredGrid>\n");
1934     /* write the xml footer */
1935     fprintf(fileHandle_p, "</VTKFile>\n");
1936     /* close the file */
1937     fclose(fileHandle_p);
1938     return;
1939     }
1940 dhawcroft 793 #endif
1941    

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26