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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26