/[escript]/branches/domexper/dudley/src/Mesh_saveVTK.c
ViewVC logotype

Annotation of /branches/domexper/dudley/src/Mesh_saveVTK.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 812 - (hide annotations)
Thu Aug 17 05:49:59 2006 UTC (12 years, 9 months ago) by dhawcroft
Original Path: trunk/finley/src/Mesh_saveVTK.c
File MIME type: text/plain
File size: 61225 byte(s)


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26