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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26