/[escript]/branches/domexper/dudley/src/Mesh_saveVTK.c
ViewVC logotype

Contents of /branches/domexper/dudley/src/Mesh_saveVTK.c

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26