/[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 812 - (show annotations)
Thu Aug 17 05:49:59 2006 UTC (13 years, 1 month ago) by dhawcroft
File MIME type: text/plain
File size: 61225 byte(s)


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26