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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26