/[escript]/trunk/finley/src/Mesh_saveVTK.c
ViewVC logotype

Contents of /trunk/finley/src/Mesh_saveVTK.c

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26