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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26