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

Contents of /trunk-mpi-branch/finley/src/Mesh_saveVTK.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 968 - (show annotations)
Tue Feb 13 22:56:57 2007 UTC (13 years, 3 months ago) by ksteube
File MIME type: text/plain
File size: 60873 byte(s)
Branch for MPI solution of implicit problems
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 if (isExpanded(data_pp[i_data])) {
870 rtmp = 0.;
871 for (j=0; j<numPointsPerSample; j++) rtmp += values[INDEX2(n,j,nComp)];
872 sampleAvg[n] = rtmp/numPointsPerSample;
873 } else {
874 sampleAvg[n] = values[n];
875 }
876 }
877 // if the number of required components is more than the number
878 // of actual components, pad with zeros
879
880 // probably only need to get shape of first element
881 // write the data different ways for scalar, vector and tensor
882 if (nCompReqd == 1)
883 {
884 sprintf(tmpbuf, " %e", sampleAvg[0]);
885 __STRCAT(largebuf,tmpbuf,tsz)
886 }
887 else if (nCompReqd == 3)
888 {
889 // write out the data
890 for (m=0; m<shape; m++)
891 {
892 sprintf(tmpbuf, " %e", sampleAvg[m]);
893 __STRCAT(largebuf,tmpbuf,tsz)
894 }
895 for (m=0; m<nCompReqd-shape; m++)
896 {
897 __STRCAT(largebuf,zero,tsz)
898 }
899 }
900 else if (nCompReqd == 9)
901 {
902 // tensor data, so have a 3x3 matrix to output as a row
903 // of 9 data points
904 count = 0;
905 for (m=0; m<shape; m++)
906 {
907 for (n=0; n<shape; n++)
908 {
909 sprintf(tmpbuf, " %e", sampleAvg[count]);
910 __STRCAT(largebuf,tmpbuf,tsz)
911 count++;
912 }
913 for (n=0; n<3-shape; n++)
914 {
915 __STRCAT(largebuf,zero,tsz)
916 }
917 }
918 for (m=0; m<3-shape; m++)
919 for (n=0; n<3; n++)
920 {
921 __STRCAT(largebuf,zero,tsz)
922 }
923 }
924 __STRCAT(largebuf,newline,tsz)
925 }
926 largebuf[tsz] = '\0';
927 MPI_File_write_ordered(fh,largebuf,tsz,MPI_CHAR,&status);
928 MEMFREE(largebuf);
929 if( myRank == 0)
930 {
931 char *tag = "</DataArray>\n";
932 MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);
933 MPI_Wait(&req,&status);
934 }
935
936 }
937 }
938 // closing celldata tag
939 if(myRank == 0)
940 {
941 char* tag = "</CellData>\n";
942 MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);
943 MPI_Wait(&req,&status);
944 }
945
946 MPIO_DEBUG(" Done Writing Cell Data ")
947 }
948
949
950 // Write Point Data Header Tags
951 if( myRank == 0)
952 {
953 char header[600];
954 char tmpBuf[50];
955
956 if (write_pointdata)
957 {
958 MPIO_DEBUG(" Writing Pointdata... ")
959 // mark the active data arrays
960 bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
961 sprintf(header, "<PointData");
962 for (i_data =0 ;i_data<num_data;++i_data)
963 {
964 if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])
965 {
966 // if the rank == 0: --> scalar data
967 // if the rank == 1: --> vector data
968 // if the rank == 2: --> tensor data
969
970 switch(getDataPointRank(data_pp[i_data]))
971 {
972 case 0:
973 if (! set_scalar)
974 {
975 sprintf(tmpBuf," Scalars=\"%s\"",names_p[i_data]);
976 strcat(header,tmpBuf);
977 set_scalar=TRUE;
978 }
979 break;
980 case 1:
981 if (! set_vector)
982 {
983 sprintf(tmpBuf," Vectors=\"%s\"",names_p[i_data]);
984 strcat(header,tmpBuf);
985 set_vector=TRUE;
986 }
987 break;
988 case 2:
989 if (! set_tensor)
990 {
991 sprintf(tmpBuf," Tensors=\"%s\"",names_p[i_data]);
992 strcat(header,tmpBuf);
993 set_tensor=TRUE;
994 }
995 break;
996 default:
997 sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);
998 Finley_setError(VALUE_ERROR,error_msg);
999 return;
1000 }
1001 }
1002 }
1003 strcat(header, ">\n");
1004 MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);
1005 MPI_Wait(&req,&status);
1006 }
1007 }
1008
1009 // write actual data
1010 if(write_pointdata)
1011 {
1012 for (i_data =0 ;i_data<num_data;++i_data)
1013 {
1014 if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])
1015 {
1016 numPointsPerSample = elements->ReferenceElement->numQuadNodes;
1017 rank = getDataPointRank(data_pp[i_data]);
1018 nComp = getDataPointSize(data_pp[i_data]);
1019 nCompReqd=1; // the number of components required by vtk
1020 shape=0;
1021 if (rank == 0)
1022 {
1023 nCompReqd = 1;
1024 }
1025 else if (rank == 1)
1026 {
1027 shape=getDataPointShape(data_pp[i_data], 0);
1028 if (shape>3)
1029 {
1030 Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");
1031 return;
1032 }
1033 nCompReqd = 3;
1034 }
1035 else
1036 {
1037 shape=getDataPointShape(data_pp[i_data], 0);
1038 if (shape>3 || shape != getDataPointShape(data_pp[i_data], 1))
1039 {
1040 Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
1041 return;
1042 }
1043 nCompReqd = 9;
1044 }
1045
1046 if( myRank == 0)
1047 {
1048 char header[250];
1049 sprintf(header, "<DataArray Name=\"%s\" type=\"Float64\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);
1050 MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req);
1051 MPI_Wait(&req,&status);
1052 }
1053 // write out the data
1054 // if the number of required components is more than the number
1055 // of actual components, pad with zeros
1056
1057 char tmpbuf[15];
1058 char* largebuf = MEMALLOC(nCompReqd*numLocalNodes*15 + numLocal*nCompReqd + nCompReqd + 15,char);
1059 largebuf[0] = '\0';
1060 bool_t do_write=TRUE;
1061 size_t tsz = 0;
1062
1063 for(k=0;k < nodeCache.size;k++)
1064 {
1065 i = nodeCache.values[k];
1066
1067 if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1068 {
1069 if (mesh_p->Nodes->toReduced[i]>=0)
1070 {
1071 switch(getFunctionSpaceType(data_pp[i_data]))
1072 {
1073 case FINLEY_DEGREES_OF_FREEDOM:
1074 values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);
1075 break;
1076 case FINLEY_REDUCED_DEGREES_OF_FREEDOM:
1077 values = getSampleData(data_pp[i_data],mesh_p->Nodes->reducedDegreeOfFreedom[i]);
1078 break;
1079 case FINLEY_NODES:
1080 values = getSampleData(data_pp[i_data],i);
1081 break;
1082 }
1083 do_write=TRUE;
1084 }
1085 else
1086 {
1087 do_write=FALSE;
1088 }
1089 }
1090 else
1091 {
1092 do_write=TRUE;
1093 switch(getFunctionSpaceType(data_pp[i_data]))
1094 {
1095 case FINLEY_DEGREES_OF_FREEDOM:
1096 values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);
1097 break;
1098 case FINLEY_NODES:
1099 values = getSampleData(data_pp[i_data],i);
1100 break;
1101 }
1102 }
1103 // write the data different ways for scalar, vector and tensor
1104 if (do_write)
1105 {
1106 if (nCompReqd == 1)
1107 {
1108 sprintf(tmpbuf," %e", values[0]);
1109 __STRCAT(largebuf,tmpbuf,tsz)
1110 }
1111 else if (nCompReqd == 3)
1112 {
1113 for (m=0; m<shape; m++)
1114 {
1115
1116 sprintf(tmpbuf," %e",values[m]);
1117 __STRCAT(largebuf,tmpbuf,tsz)
1118 }
1119 for (m=0; m<nCompReqd-shape; m++)
1120 {
1121 __STRCAT(largebuf,zero,tsz)
1122 }
1123 }
1124 else if (nCompReqd == 9)
1125 {
1126 // tensor data, so have a 3x3 matrix to output as a row
1127 // of 9 data points
1128 count = 0;
1129 for (m=0; m<shape; m++)
1130 {
1131 for (n=0; n<shape; n++)
1132 {
1133 sprintf(tmpbuf," %e", values[count]);
1134 __STRCAT(largebuf,tmpbuf,tsz)
1135 count++;
1136 }
1137 for (n=0; n<3-shape; n++)
1138 {
1139 __STRCAT(largebuf,zero,tsz)
1140 }
1141 }
1142 for (m=0; m<3-shape; m++)
1143 {
1144 for (n=0; n<3; n++)
1145 {
1146 __STRCAT(largebuf,zero,tsz)
1147 }
1148 }
1149 }
1150 __STRCAT(largebuf,newline,tsz)
1151 }
1152
1153 }
1154 // Write out local data
1155
1156 largebuf[tsz] = '\0';
1157 MPI_File_write_ordered(fh,largebuf,tsz,MPI_CHAR,&status);
1158 MEMFREE(largebuf);
1159 if( myRank == 0)
1160 {
1161 char *tag = "</DataArray>\n";
1162 MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);
1163 MPI_Wait(&req,&status);
1164 }
1165 }
1166 }
1167 // Finish off with closing tag
1168 if(myRank == 0)
1169 {
1170 char* tag = "</PointData>\n";
1171 MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req);
1172 MPI_Wait(&req,&status);
1173 }
1174 MPIO_DEBUG(" Done Writing Pointdata ")
1175 }
1176 // end write_pointdata
1177
1178 // tag and bag...
1179 if (myRank == 0)
1180 {
1181 char *footer = "</Piece>\n</UnstructuredGrid>\n</VTKFile>";
1182 MPI_File_iwrite_shared(fh,footer,strlen(footer),MPI_CHAR,&req);
1183 MPI_Wait(&req,&status);
1184 }
1185
1186 MEMFREE(nodesGlobal);
1187 MEMFREE(nodeCache.values);
1188 MEMFREE(elementCache.values);
1189 #ifdef MPIO_HINTS
1190 MPI_Info_free(&infoHints);
1191 #undef MPIO_HINTS
1192 #endif
1193 MPI_File_close(&fh);
1194 MPIO_DEBUG(" ***** Exit saveVTK ***** ")
1195 #undef __STRCAT
1196 }
1197
1198 #undef MPIO_DEBUG
1199 #else
1200
1201
1202
1203
1204 void Finley_Mesh_saveVTK(const char * filename_p, Finley_Mesh *mesh_p, const dim_t num_data,char* *names_p, escriptDataC* *data_pp)
1205 {
1206 char error_msg[LenErrorMsg_MAX];
1207 /* if there is no mesh we just return */
1208 if (mesh_p==NULL) return;
1209
1210 int i, j, k, numVTKNodesPerElement,i_data,m, count, n, rank,shape, numPoints, cellType, numCells,
1211 nDim, numPointsPerSample, nComp, nCompReqd;
1212
1213 index_t j2;
1214 double* values, rtmp;
1215 char elemTypeStr[32];
1216
1217 /* open the file and check handle */
1218
1219 FILE * fileHandle_p = fopen(filename_p, "w");
1220 if (fileHandle_p==NULL)
1221 {
1222 sprintf(error_msg, "saveVTK: File %s could not be opened for writing.", filename_p);
1223 Finley_setError(IO_ERROR,error_msg);
1224 return;
1225 }
1226 /* find the mesh type to be written */
1227 int nodetype=FINLEY_DEGREES_OF_FREEDOM;
1228 int elementtype=FINLEY_UNKNOWN;
1229 bool_t isCellCentered[num_data],write_celldata=FALSE,write_pointdata=FALSE;
1230 for (i_data=0;i_data<num_data;++i_data)
1231 {
1232 if (! isEmpty(data_pp[i_data]))
1233 {
1234 switch(getFunctionSpaceType(data_pp[i_data]))
1235 {
1236 case FINLEY_DEGREES_OF_FREEDOM:
1237 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1238 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)
1239 {
1240 elementtype=FINLEY_ELEMENTS;
1241 }
1242 else
1243 {
1244 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1245 fclose(fileHandle_p);
1246 return;
1247 }
1248 isCellCentered[i_data]=FALSE;
1249 break;
1250 case FINLEY_REDUCED_DEGREES_OF_FREEDOM:
1251 nodetype = FINLEY_REDUCED_DEGREES_OF_FREEDOM;
1252 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)
1253 {
1254 elementtype=FINLEY_ELEMENTS;
1255 }
1256 else
1257 {
1258 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1259 fclose(fileHandle_p);
1260 return;
1261 }
1262 isCellCentered[i_data]=FALSE;
1263 break;
1264 case FINLEY_NODES:
1265 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1266 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)
1267 {
1268 elementtype=FINLEY_ELEMENTS;
1269 }
1270 else
1271 {
1272 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1273 fclose(fileHandle_p);
1274 return;
1275 }
1276 isCellCentered[i_data]=FALSE;
1277 break;
1278 case FINLEY_ELEMENTS:
1279 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1280 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS)
1281 {
1282 elementtype=FINLEY_ELEMENTS;
1283 }
1284 else
1285 {
1286 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1287 fclose(fileHandle_p);
1288 return;
1289 }
1290 isCellCentered[i_data]=TRUE;
1291 break;
1292 case FINLEY_FACE_ELEMENTS:
1293 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1294 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_FACE_ELEMENTS)
1295 {
1296 elementtype=FINLEY_FACE_ELEMENTS;
1297 }
1298 else
1299 {
1300 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1301 fclose(fileHandle_p);
1302 return;
1303 }
1304 isCellCentered[i_data]=TRUE;
1305 break;
1306 case FINLEY_POINTS:
1307 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1308 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_POINTS)
1309 {
1310 elementtype=FINLEY_POINTS;
1311 }
1312 else
1313 {
1314 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1315 fclose(fileHandle_p);
1316 return;
1317 }
1318 isCellCentered[i_data]=TRUE;
1319 break;
1320 case FINLEY_CONTACT_ELEMENTS_1:
1321 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1322 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1)
1323 {
1324 elementtype=FINLEY_CONTACT_ELEMENTS_1;
1325 }
1326 else
1327 {
1328 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1329 fclose(fileHandle_p);
1330 return;
1331 }
1332 isCellCentered[i_data]=TRUE;
1333 break;
1334 case FINLEY_CONTACT_ELEMENTS_2:
1335 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
1336 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1)
1337 {
1338 elementtype=FINLEY_CONTACT_ELEMENTS_1;
1339 }
1340 else
1341 {
1342 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
1343 fclose(fileHandle_p);
1344 return;
1345 }
1346 isCellCentered[i_data]=TRUE;
1347 break;
1348 default:
1349 sprintf(error_msg,"saveVTK: Finley does not know anything about function space type %d",getFunctionSpaceType(data_pp[i_data]));
1350 Finley_setError(TYPE_ERROR,error_msg);
1351 fclose(fileHandle_p);
1352 return;
1353 }
1354 if (isCellCentered[i_data])
1355 {
1356 write_celldata=TRUE;
1357 }
1358 else
1359 {
1360 write_pointdata=TRUE;
1361 }
1362 }
1363 }
1364 /* select nomber of points and the mesh component */
1365 numPoints = mesh_p->Nodes->numNodes;
1366 if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1367 {
1368 numPoints = mesh_p->Nodes->reducedNumNodes;
1369 }
1370 else
1371 {
1372 numPoints = mesh_p->Nodes->numNodes;
1373 }
1374 if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;
1375 Finley_ElementFile* elements=NULL;
1376 switch(elementtype)
1377 {
1378 case FINLEY_ELEMENTS:
1379 elements=mesh_p->Elements;
1380 break;
1381 case FINLEY_FACE_ELEMENTS:
1382 elements=mesh_p->FaceElements;
1383 break;
1384 case FINLEY_POINTS:
1385 elements=mesh_p->Points;
1386 break;
1387 case FINLEY_CONTACT_ELEMENTS_1:
1388 elements=mesh_p->ContactElements;
1389 break;
1390 }
1391 if (elements==NULL)
1392 {
1393 Finley_setError(SYSTEM_ERROR,"saveVTK: undefined element file");
1394 fclose(fileHandle_p);
1395 return;
1396 }
1397 /* map finley element type to VTK element type */
1398 numCells = elements->numElements;
1399 ElementTypeId TypeId;
1400 if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1401 {
1402 TypeId = elements->LinearReferenceElement->Type->TypeId;
1403 }
1404 else
1405 {
1406 TypeId = elements->ReferenceElement->Type->TypeId;
1407 }
1408
1409 switch(TypeId)
1410 {
1411 case Point1:
1412 case Line2Face:
1413 case Line3Face:
1414 case Point1_Contact:
1415 case Line2Face_Contact:
1416 case Line3Face_Contact:
1417 cellType = VTK_VERTEX;
1418 numVTKNodesPerElement = 1;
1419 strcpy(elemTypeStr, "VTK_VERTEX");
1420 break;
1421
1422 case Line2:
1423 case Tri3Face:
1424 case Rec4Face:
1425 case Line2_Contact:
1426 case Tri3_Contact:
1427 case Tri3Face_Contact:
1428 case Rec4Face_Contact:
1429 cellType = VTK_LINE;
1430 numVTKNodesPerElement = 2;
1431 strcpy(elemTypeStr, "VTK_LINE");
1432 break;
1433
1434 case Tri3:
1435 case Tet4Face:
1436 case Tet4Face_Contact:
1437 cellType = VTK_TRIANGLE;
1438 numVTKNodesPerElement = 3;
1439 strcpy(elemTypeStr, "VTK_TRIANGLE");
1440 break;
1441
1442 case Rec4:
1443 case Hex8Face:
1444 case Rec4_Contact:
1445 case Hex8Face_Contact:
1446 cellType = VTK_QUAD;
1447 numVTKNodesPerElement = 4;
1448 strcpy(elemTypeStr, "VTK_QUAD");
1449 break;
1450
1451 case Tet4:
1452 cellType = VTK_TETRA;
1453 numVTKNodesPerElement = 4;
1454 strcpy(elemTypeStr, "VTK_TETRA");
1455 break;
1456
1457 case Hex8:
1458 cellType = VTK_HEXAHEDRON;
1459 numVTKNodesPerElement = 8;
1460 strcpy(elemTypeStr, "VTK_HEXAHEDRON");
1461 break;
1462
1463 case Line3:
1464 case Tri6Face:
1465 case Rec8Face:
1466 case Line3_Contact:
1467 case Tri6Face_Contact:
1468 case Rec8Face_Contact:
1469 cellType = VTK_QUADRATIC_EDGE;
1470 numVTKNodesPerElement = 3;
1471 strcpy(elemTypeStr, "VTK_QUADRATIC_EDGE");
1472 break;
1473
1474 case Tri6:
1475 case Tet10Face:
1476 case Tri6_Contact:
1477 case Tet10Face_Contact:
1478 cellType = VTK_QUADRATIC_TRIANGLE;
1479 numVTKNodesPerElement = 6;
1480 strcpy(elemTypeStr, "VTK_QUADRATIC_TRIANGLE");
1481 break;
1482
1483 case Rec8:
1484 case Hex20Face:
1485 case Rec8_Contact:
1486 case Hex20Face_Contact:
1487 cellType = VTK_QUADRATIC_QUAD;
1488 numVTKNodesPerElement = 8;
1489 strcpy(elemTypeStr, "VTK_QUADRATIC_QUAD");
1490 break;
1491
1492 case Tet10:
1493 cellType = VTK_QUADRATIC_TETRA;
1494 numVTKNodesPerElement = 10;
1495 strcpy(elemTypeStr, "VTK_QUADRATIC_TETRA");
1496 break;
1497
1498 case Hex20:
1499 cellType = VTK_QUADRATIC_HEXAHEDRON;
1500 numVTKNodesPerElement = 20;
1501 strcpy(elemTypeStr, "VTK_QUADRATIC_HEXAHEDRON");
1502 break;
1503
1504 default:
1505 sprintf(error_msg, "saveVTK: Element type %s is not supported by VTK",elements->ReferenceElement->Type->Name);
1506 Finley_setError(VALUE_ERROR,error_msg);
1507 fclose(fileHandle_p);
1508 return;
1509 }
1510 /* xml header */
1511 fprintf(fileHandle_p, "<?xml version=\"1.0\"?>\n");
1512 fprintf(fileHandle_p, "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\">\n");
1513
1514 /* finley uses an unstructured mesh, so UnstructuredGrid *should* work */
1515 fprintf(fileHandle_p, "<UnstructuredGrid>\n");
1516
1517 /* is there only one "piece" to the data?? */
1518 fprintf(fileHandle_p, "<Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n",numPoints, numCells);
1519 /* now for the points; equivalent to positions section in saveDX() */
1520 /* "The points element explicitly defines coordinates for each point
1521 * individually. It contains one DataArray element describing an array
1522 * with three components per value, each specifying the coordinates of one
1523 * point" - from Vtk User's Guide
1524 */
1525 fprintf(fileHandle_p, "<Points>\n");
1526 /*
1527 * the reason for this if statement is explained in the long comment below
1528 */
1529 nDim = mesh_p->Nodes->numDim;
1530 fprintf(fileHandle_p, "<DataArray NumberOfComponents=\"%d\" type=\"Float64\" format=\"ascii\">\n",MAX(3,nDim));
1531 /* vtk/mayavi doesn't like 2D data, it likes 3D data with a degenerate
1532 * third dimension to handle 2D data (like a sheet of paper). So, if
1533 * nDim is 2, we have to append zeros to the array to get this third
1534 * dimension, and keep the visualisers happy.
1535 * Indeed, if nDim is less than 3, must pad all empty dimensions, so
1536 * that the total number of dims is 3.
1537 */
1538 if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1539 {
1540 for (i = 0; i < mesh_p->Nodes->numNodes; i++)
1541 {
1542 if (mesh_p->Nodes->toReduced[i]>=0)
1543 {
1544 for (j = 0; j < nDim; j++) fprintf(fileHandle_p, " %e",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);
1545 for (k=0; k<3-nDim; k++) fprintf(fileHandle_p, " %e",0.);
1546 fprintf(fileHandle_p, "\n");
1547 }
1548 }
1549 }
1550 else
1551 {
1552 for (i = 0; i < mesh_p->Nodes->numNodes; i++)
1553 {
1554
1555 for (j = 0; j < nDim; j++) fprintf(fileHandle_p, " %e",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);
1556 for (k=0; k<3-nDim; k++) fprintf(fileHandle_p, " %e",0.);
1557 fprintf(fileHandle_p, "\n");
1558 }
1559
1560 }
1561 fprintf(fileHandle_p, "</DataArray>\n");
1562 fprintf(fileHandle_p, "</Points>\n");
1563
1564 /* write out the DataArray element for the connectivity */
1565
1566 int NN = elements->ReferenceElement->Type->numNodes;
1567 fprintf(fileHandle_p, "<Cells>\n");
1568 fprintf(fileHandle_p, "<DataArray Name=\"connectivity\" type=\"Int32\" format=\"ascii\">\n");
1569
1570 if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1571 {
1572 for (i = 0; i < numCells; i++)
1573 {
1574 for (j = 0; j < numVTKNodesPerElement; j++)
1575 fprintf(fileHandle_p,"%d ",mesh_p->Nodes->toReduced[elements->Nodes[INDEX2(elements->ReferenceElement->Type->linearNodes[j], i, NN)]]);
1576 fprintf(fileHandle_p, "\n");
1577 }
1578 }
1579 else if (VTK_QUADRATIC_HEXAHEDRON==cellType)
1580 {
1581 for (i = 0; i < numCells; i++)
1582 {
1583 fprintf(fileHandle_p,"%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d\n",
1584 elements->Nodes[INDEX2(0, i, NN)],
1585 elements->Nodes[INDEX2(1, i, NN)],
1586 elements->Nodes[INDEX2(2, i, NN)],
1587 elements->Nodes[INDEX2(3, i, NN)],
1588 elements->Nodes[INDEX2(4, i, NN)],
1589 elements->Nodes[INDEX2(5, i, NN)],
1590 elements->Nodes[INDEX2(6, i, NN)],
1591 elements->Nodes[INDEX2(7, i, NN)],
1592 elements->Nodes[INDEX2(8, i, NN)],
1593 elements->Nodes[INDEX2(9, i, NN)],
1594 elements->Nodes[INDEX2(10, i, NN)],
1595 elements->Nodes[INDEX2(11, i, NN)],
1596 elements->Nodes[INDEX2(16, i, NN)],
1597 elements->Nodes[INDEX2(17, i, NN)],
1598 elements->Nodes[INDEX2(18, i, NN)],
1599 elements->Nodes[INDEX2(19, i, NN)],
1600 elements->Nodes[INDEX2(12, i, NN)],
1601 elements->Nodes[INDEX2(13, i, NN)],
1602 elements->Nodes[INDEX2(14, i, NN)],
1603 elements->Nodes[INDEX2(15, i, NN)]);
1604 }
1605 }
1606 else if (numVTKNodesPerElement!=NN)
1607 {
1608 for (i = 0; i < numCells; i++)
1609 {
1610 for (j = 0; j < numVTKNodesPerElement; j++) fprintf(fileHandle_p,"%d ", elements->Nodes[INDEX2(elements->ReferenceElement->Type->geoNodes[j], i, NN)]);
1611 fprintf(fileHandle_p, "\n");
1612 }
1613 }
1614 else
1615 {
1616 for (i = 0; i < numCells; i++)
1617 {
1618 for (j = 0; j < numVTKNodesPerElement; j++) fprintf(fileHandle_p,"%d ", elements->Nodes[INDEX2(j, i, NN)]);
1619 fprintf(fileHandle_p, "\n");
1620 }
1621 }
1622 fprintf(fileHandle_p, "</DataArray>\n");
1623
1624 /* write out the DataArray element for the offsets */
1625 fprintf(fileHandle_p, "<DataArray Name=\"offsets\" type=\"Int32\" format=\"ascii\">\n");
1626 for (i=numVTKNodesPerElement; i<=numCells*numVTKNodesPerElement; i+=numVTKNodesPerElement) fprintf(fileHandle_p, "%d\n", i);
1627 fprintf(fileHandle_p, "</DataArray>\n");
1628
1629 /* write out the DataArray element for the types */
1630 fprintf(fileHandle_p, "<DataArray Name=\"types\" type=\"UInt8\" format=\"ascii\">\n");
1631 for (i=0; i<numCells; i++) fprintf(fileHandle_p, "%d\n", cellType);
1632 fprintf(fileHandle_p, "</DataArray>\n");
1633
1634 /* finish off the <Cells> element */
1635 fprintf(fileHandle_p, "</Cells>\n");
1636
1637 /* cell data */
1638 if (write_celldata)
1639 {
1640 /* mark the active data arrays */
1641 bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
1642 fprintf(fileHandle_p, "<CellData");
1643 for (i_data =0 ;i_data<num_data;++i_data)
1644 {
1645 if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data])
1646 {
1647 /* if the rank == 0: --> scalar data
1648 * if the rank == 1: --> vector data
1649 * if the rank == 2: --> tensor data
1650 */
1651 switch(getDataPointRank(data_pp[i_data]))
1652 {
1653 case 0:
1654 if (! set_scalar)
1655 {
1656 fprintf(fileHandle_p," Scalars=\"%s\"",names_p[i_data]);
1657 set_scalar=TRUE;
1658 }
1659 break;
1660 case 1:
1661 if (! set_vector)
1662 {
1663 fprintf(fileHandle_p," Vectors=\"%s\"",names_p[i_data]);
1664 set_vector=TRUE;
1665 }
1666 break;
1667 case 2:
1668 if (! set_tensor)
1669 {
1670 fprintf(fileHandle_p," Tensors=\"%s\"",names_p[i_data]);
1671 set_tensor=TRUE;
1672 }
1673 break;
1674 default:
1675 sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);
1676 Finley_setError(VALUE_ERROR,error_msg);
1677 fclose(fileHandle_p);
1678 return;
1679 }
1680 }
1681 }
1682 fprintf(fileHandle_p, ">\n");
1683 /* write the arrays */
1684 for (i_data =0 ;i_data<num_data;++i_data)
1685 {
1686 if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data])
1687 {
1688 numPointsPerSample = elements->ReferenceElement->numQuadNodes;
1689 rank = getDataPointRank(data_pp[i_data]);
1690 nComp = getDataPointSize(data_pp[i_data]);
1691 nCompReqd=1; /* the number of components required by vtk */
1692 shape=0;
1693 if (rank == 0)
1694 {
1695 nCompReqd = 1;
1696 }
1697 else if (rank == 1)
1698 {
1699 shape=getDataPointShape(data_pp[i_data], 0);
1700 if (shape>3)
1701 {
1702 Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");
1703 fclose(fileHandle_p);
1704 return;
1705 }
1706 nCompReqd = 3;
1707 }
1708 else
1709 {
1710 shape=getDataPointShape(data_pp[i_data], 0);
1711 if (shape>3 || shape != getDataPointShape(data_pp[i_data], 1))
1712 {
1713 Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
1714 fclose(fileHandle_p);
1715 return;
1716 }
1717 nCompReqd = 9;
1718 }
1719 fprintf(fileHandle_p, "<DataArray Name=\"%s\" type=\"Float64\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);
1720
1721 double sampleAvg[nComp];
1722 for (i=0; i<numCells; i++)
1723 {
1724 values = getSampleData(data_pp[i_data], i);
1725 /* averaging over the number of points in the sample */
1726 for (k=0; k<nComp; k++)
1727 {
1728 if (isExpanded(data_pp[i_data])) {
1729 rtmp = 0.;
1730 for (j=0; j<numPointsPerSample; j++) rtmp += values[INDEX2(k,j,nComp)];
1731 sampleAvg[k] = rtmp/numPointsPerSample;
1732 } else {
1733 sampleAvg[k] = values[k];
1734 }
1735
1736 }
1737 /* if the number of required components is more than the number
1738 * of actual components, pad with zeros
1739 */
1740 /* probably only need to get shape of first element */
1741 /* write the data different ways for scalar, vector and tensor */
1742 if (nCompReqd == 1)
1743 {
1744 fprintf(fileHandle_p, " %e", sampleAvg[0]);
1745 }
1746 else if (nCompReqd == 3)
1747 {
1748 /* write out the data */
1749 for (m=0; m<shape; m++) fprintf(fileHandle_p, " %e", sampleAvg[m]);
1750 for (m=0; m<nCompReqd-shape; m++) fprintf(fileHandle_p, " %e", 0.);
1751 }
1752 else if (nCompReqd == 9)
1753 {
1754 /* tensor data, so have a 3x3 matrix to output as a row
1755 * of 9 data points */
1756 count = 0;
1757 for (m=0; m<shape; m++)
1758 {
1759 for (n=0; n<shape; n++)
1760 {
1761 fprintf(fileHandle_p, " %e", sampleAvg[count]);
1762 count++;
1763 }
1764 for (n=0; n<3-shape; n++) fprintf(fileHandle_p, " %e", 0.);
1765 }
1766 for (m=0; m<3-shape; m++)
1767 for (n=0; n<3; n++) fprintf(fileHandle_p, " %e", 0.);
1768 }
1769 fprintf(fileHandle_p, "\n");
1770 }
1771 fprintf(fileHandle_p, "</DataArray>\n");
1772 }
1773 }
1774 fprintf(fileHandle_p, "</CellData>\n");
1775 }
1776 /* point data */
1777 if (write_pointdata)
1778 {
1779 /* mark the active data arrays */
1780 bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
1781 fprintf(fileHandle_p, "<PointData");
1782 for (i_data =0 ;i_data<num_data;++i_data)
1783 {
1784 if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])
1785 {
1786 /* if the rank == 0: --> scalar data
1787 * if the rank == 1: --> vector data
1788 * if the rank == 2: --> tensor data
1789 */
1790 switch(getDataPointRank(data_pp[i_data]))
1791 {
1792 case 0:
1793 if (! set_scalar)
1794 {
1795 fprintf(fileHandle_p," Scalars=\"%s\"",names_p[i_data]);
1796 set_scalar=TRUE;
1797 }
1798 break;
1799 case 1:
1800 if (! set_vector)
1801 {
1802 fprintf(fileHandle_p," Vectors=\"%s\"",names_p[i_data]);
1803 set_vector=TRUE;
1804 }
1805 break;
1806 case 2:
1807 if (! set_tensor)
1808 {
1809 fprintf(fileHandle_p," Tensors=\"%s\"",names_p[i_data]);
1810 set_tensor=TRUE;
1811 }
1812 break;
1813 default:
1814 sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);
1815 Finley_setError(VALUE_ERROR,error_msg);
1816 fclose(fileHandle_p);
1817 return;
1818 }
1819 }
1820 }
1821 fprintf(fileHandle_p, ">\n");
1822 /* write the arrays */
1823 for (i_data =0 ;i_data<num_data;++i_data)
1824 {
1825 if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data])
1826 {
1827 numPointsPerSample = elements->ReferenceElement->numQuadNodes;
1828 rank = getDataPointRank(data_pp[i_data]);
1829 nComp = getDataPointSize(data_pp[i_data]);
1830 nCompReqd=1; /* the number of components required by vtk */
1831 shape=0;
1832 if (rank == 0)
1833 {
1834 nCompReqd = 1;
1835 }
1836 else if (rank == 1)
1837 {
1838 shape=getDataPointShape(data_pp[i_data], 0);
1839 if (shape>3)
1840 {
1841 Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");
1842 fclose(fileHandle_p);
1843 return;
1844 }
1845 nCompReqd = 3;
1846 }
1847 else
1848 {
1849 shape=getDataPointShape(data_pp[i_data], 0);
1850 if (shape>3 || shape != getDataPointShape(data_pp[i_data], 1))
1851 {
1852 Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
1853 fclose(fileHandle_p);
1854 return;
1855 }
1856 nCompReqd = 9;
1857 }
1858 fprintf(fileHandle_p, "<DataArray Name=\"%s\" type=\"Float64\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);
1859 /* write out the data */
1860 /* if the number of required components is more than the number
1861 * of actual components, pad with zeros
1862 */
1863 bool_t do_write=TRUE;
1864 for (i=0; i<mesh_p->Nodes->numNodes; i++)
1865 {
1866 if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM)
1867 {
1868 if (mesh_p->Nodes->toReduced[i]>=0)
1869 {
1870 switch(getFunctionSpaceType(data_pp[i_data]))
1871 {
1872 case FINLEY_DEGREES_OF_FREEDOM:
1873 values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);
1874 break;
1875 case FINLEY_REDUCED_DEGREES_OF_FREEDOM:
1876 values = getSampleData(data_pp[i_data],mesh_p->Nodes->reducedDegreeOfFreedom[i]);
1877 break;
1878 case FINLEY_NODES:
1879 values = getSampleData(data_pp[i_data],i);
1880 break;
1881 }
1882 do_write=TRUE;
1883 }
1884 else
1885 {
1886 do_write=FALSE;
1887 }
1888 }
1889 else
1890 {
1891 do_write=TRUE;
1892 switch(getFunctionSpaceType(data_pp[i_data]))
1893 {
1894 case FINLEY_DEGREES_OF_FREEDOM:
1895 values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);
1896 break;
1897 case FINLEY_NODES:
1898 values = getSampleData(data_pp[i_data],i);
1899 break;
1900 }
1901 }
1902 /* write the data different ways for scalar, vector and tensor */
1903 if (do_write)
1904 {
1905 if (nCompReqd == 1)
1906 {
1907 fprintf(fileHandle_p, " %e", values[0]);
1908 }
1909 else if (nCompReqd == 3)
1910 {
1911 for (m=0; m<shape; m++) fprintf(fileHandle_p, " %e", values[m]);
1912 for (m=0; m<nCompReqd-shape; m++) fprintf(fileHandle_p, " %e", 0.);
1913 }
1914 else if (nCompReqd == 9)
1915 {
1916 /* tensor data, so have a 3x3 matrix to output as a row
1917 * of 9 data points */
1918 count = 0;
1919 for (m=0; m<shape; m++)
1920 {
1921 for (n=0; n<shape; n++)
1922 {
1923 fprintf(fileHandle_p, " %e", values[count]);
1924 count++;
1925 }
1926 for (n=0; n<3-shape; n++) fprintf(fileHandle_p, " %e", 0.);
1927 }
1928 for (m=0; m<3-shape; m++)
1929 for (n=0; n<3; n++) fprintf(fileHandle_p, " %e", 0.);
1930 }
1931 fprintf(fileHandle_p, "\n");
1932 }
1933 }
1934 fprintf(fileHandle_p, "</DataArray>\n");
1935 }
1936 }
1937 fprintf(fileHandle_p, "</PointData>\n");
1938 }
1939 /* finish off the piece */
1940 fprintf(fileHandle_p, "</Piece>\n");
1941
1942 fprintf(fileHandle_p, "</UnstructuredGrid>\n");
1943 /* write the xml footer */
1944 fprintf(fileHandle_p, "</VTKFile>\n");
1945 /* close the file */
1946 fclose(fileHandle_p);
1947 return;
1948 }
1949 #endif
1950

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26