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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26