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

Contents of /trunk/finley/src/finley/Mesh_saveVTK.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 150 - (show annotations)
Thu Sep 15 03:44:45 2005 UTC (13 years, 11 months ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Mesh_saveVTK.c
File MIME type: text/plain
File size: 23534 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-15

1 /*
2 ******************************************************************************
3 * *
4 * COPYRIGHT ACcESS 2003,2004,2005 - All Rights Reserved *
5 * *
6 * This software is the property of ACcESS. No part of this code *
7 * may be copied in any form or by any means without the expressed written *
8 * consent of ACcESS. Copying, use or modification of this software *
9 * by any unauthorised person is illegal unless that person has a software *
10 * license agreement with ACcESS. *
11 * *
12 ******************************************************************************
13 */
14
15
16 /**************************************************************/
17
18 /* writes data and mesh in a vtk file */
19
20 /**************************************************************/
21
22 /* Author: Paul Cochrane, cochrane@esscc.uq.edu.au */
23 /* Version: $Id$ */
24
25 /**************************************************************/
26
27 #include "Mesh.h"
28 #include "vtkCellType.h" /* copied from vtk source directory !!! */
29
30 /**************************************************************/
31
32 void Finley_Mesh_saveVTK(const char * filename_p, Finley_Mesh *mesh_p, escriptDataC* data_p) {
33 char error_msg[LenErrorMsg_MAX];
34 /* if there is no mesh we just return */
35 if (mesh_p==NULL) return;
36 Finley_ElementFile* elements=NULL;
37 char elemTypeStr[32];
38 int i, j, k, numVTKNodesPerElement, isCellCentered=FALSE, nodetype=FINLEY_DEGREES_OF_FREEDOM;
39 index_t j2;
40 double* values, rtmp;
41 int nDim = mesh_p->Nodes->numDim;
42 int numPoints=0;
43
44 if (nDim==1) {
45 Finley_setError(VALUE_ERROR,"saveVTK: 1-dimensional domains are not supported.");
46 return;
47 }
48 /* get a pointer to the relevant mesh component */
49 if (isEmpty(data_p)) {
50 numPoints = mesh_p->Nodes->numNodes;
51 elements=mesh_p->Elements;
52 } else {
53 switch(getFunctionSpaceType(data_p)) {
54 case(FINLEY_DEGREES_OF_FREEDOM):
55 numPoints = mesh_p->Nodes->numNodes;
56 nodetype = FINLEY_DEGREES_OF_FREEDOM;
57 isCellCentered = FALSE;
58 elements = mesh_p->Elements;
59 break;
60 case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):
61 numPoints = mesh_p->Nodes->reducedNumNodes;
62 nodetype =FINLEY_REDUCED_DEGREES_OF_FREEDOM;
63 isCellCentered = FALSE;
64 elements = mesh_p->Elements;
65 break;
66 case(FINLEY_NODES):
67 numPoints = mesh_p->Nodes->numNodes;
68 nodetype=FINLEY_NODES;
69 isCellCentered=FALSE;
70 elements=mesh_p->Elements;
71 break;
72 case(FINLEY_ELEMENTS):
73 numPoints = mesh_p->Nodes->numNodes;
74 nodetype=FINLEY_NODES;
75 isCellCentered=TRUE;
76 elements=mesh_p->Elements;
77 break;
78 case(FINLEY_FACE_ELEMENTS):
79 numPoints = mesh_p->Nodes->numNodes;
80 nodetype=FINLEY_NODES;
81 isCellCentered=TRUE;
82 elements=mesh_p->FaceElements;
83 break;
84 case(FINLEY_POINTS):
85 numPoints = mesh_p->Nodes->numNodes;
86 nodetype=FINLEY_NODES;
87 isCellCentered=TRUE;
88 elements=mesh_p->Points;
89 break;
90 case(FINLEY_CONTACT_ELEMENTS_1):
91 case(FINLEY_CONTACT_ELEMENTS_2):
92 numPoints = mesh_p->Nodes->numNodes;
93 nodetype=FINLEY_NODES;
94 isCellCentered=TRUE;
95 elements=mesh_p->ContactElements;
96 break;
97 default:
98 sprintf(error_msg,"saveVTK: Finley does not know anything about function space type %d",getFunctionSpaceType(data_p));
99 Finley_setError(TYPE_ERROR,error_msg);
100 return;
101 }
102 }
103
104 /* the number of points */
105
106 /* the number of cells */
107 if (elements == NULL) {
108 Finley_setError(VALUE_ERROR,"saveVTK: elements object is NULL; cannot proceed");
109 return;
110 }
111 int numCells = elements->numElements;
112 /* open the file and check handle */
113 FILE * fileHandle_p = fopen(filename_p, "w");
114 if (fileHandle_p==NULL) {
115 sprintf(error_msg, "saveVTK: File %s could not be opened for writing.", filename_p);
116 Finley_setError(IO_ERROR,error_msg);
117 return;
118 }
119 /* xml header */
120 fprintf(fileHandle_p, "<?xml version=\"1.0\"?>\n");
121 fprintf(fileHandle_p, "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\">\n");
122
123 /* finley uses an unstructured mesh, so UnstructuredGrid *should* work */
124 fprintf(fileHandle_p, "<UnstructuredGrid>\n");
125
126 /* is there only one "piece" to the data?? */
127 fprintf(fileHandle_p, "<Piece "
128 "NumberOfPoints=\"%d\" "
129 "NumberOfCells=\"%d\">\n",
130 numPoints, numCells);
131 /* now for the points; equivalent to positions section in saveDX() */
132 /* "The points element explicitly defines coordinates for each point
133 * individually. It contains one DataArray element describing an array
134 * with three components per value, each specifying the coordinates of one
135 * point" - from Vtk User's Guide
136 */
137 fprintf(fileHandle_p, "<Points>\n");
138 /*
139 * the reason for this if statement is explained in the long comment below
140 */
141 if (nDim < 3) {
142 fprintf(fileHandle_p, "<DataArray "
143 "NumberOfComponents=\"3\" "
144 "type=\"Float32\" "
145 "format=\"ascii\">\n");
146 } else {
147 fprintf(fileHandle_p, "<DataArray "
148 "NumberOfComponents=\"%d\" "
149 "type=\"Float32\" "
150 "format=\"ascii\">\n",
151 nDim);
152 }
153 /* vtk/mayavi doesn't like 2D data, it likes 3D data with a degenerate
154 * third dimension to handle 2D data (like a sheet of paper). So, if
155 * nDim is 2, we have to append zeros to the array to get this third
156 * dimension, and keep the visualisers happy.
157 * Indeed, if nDim is less than 3, must pad all empty dimensions, so
158 * that the total number of dims is 3.
159 */
160 if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
161 for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
162 if (mesh_p->Nodes->toReduced[i]>=0) {
163 for (j = 0; j < nDim; j++)
164 fprintf(fileHandle_p, " %e",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);
165 for (k=0; k<3-nDim; k++) fprintf(fileHandle_p, " %e",0.);
166 fprintf(fileHandle_p, "\n");
167 }
168 }
169 } else {
170 for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
171 for (j = 0; j < nDim; j++)
172 fprintf(fileHandle_p, " %e",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);
173 for (k=0; k<3-nDim; k++) fprintf(fileHandle_p, " %e",0.);
174 fprintf(fileHandle_p, "\n");
175 }
176 }
177 fprintf(fileHandle_p, "</DataArray>\n");
178 fprintf(fileHandle_p, "</Points>\n");
179
180 /* connections */
181 /* now for the cells */
182 /* "The Cells element defines cells explicitly by specifying point
183 * connectivity and cell types. It contains three DataArray elements. The
184 * first array specifies the point connectivity. All cells' point lists
185 * are concatenated together. The second array specifies th eoffset into
186 * the connectivity array for the end of each cell. The third array
187 * specifies the type of each cell.
188 */
189 /* if no element table is present jump over the connection table */
190 if (elements!=NULL) {
191 int cellType;
192 ElementTypeId TypeId;
193
194 if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
195 TypeId = elements->LinearReferenceElement->Type->TypeId;
196 } else {
197 TypeId = elements->ReferenceElement->Type->TypeId;
198 }
199
200 switch(TypeId) {
201 case Point1:
202 cellType = VTK_VERTEX;
203 break;
204 case Line2:
205 cellType = VTK_LINE;
206 break;
207 case Line3:
208 cellType = VTK_QUADRATIC_EDGE;
209 break;
210 case Tri3:
211 cellType = VTK_TRIANGLE;
212 break;
213 case Tri6:
214 cellType = VTK_QUADRATIC_TRIANGLE;
215 break;
216 case Rec4:
217 cellType = VTK_QUAD;
218 break;
219 case Rec8:
220 cellType = VTK_QUADRATIC_QUAD;
221 break;
222 case Tet4:
223 cellType = VTK_TETRA;
224 break;
225 case Tet10:
226 cellType = VTK_QUADRATIC_TETRA;
227 break;
228 case Hex8:
229 cellType = VTK_HEXAHEDRON;
230 break;
231 case Hex20:
232 cellType = VTK_QUADRATIC_HEXAHEDRON;
233 break;
234 case Line2Face:
235 cellType = VTK_VERTEX;
236 break;
237 case Line3Face:
238 cellType = VTK_VERTEX;
239 break;
240 case Tri3Face:
241 cellType = VTK_LINE;
242 break;
243 case Tri6Face:
244 cellType = VTK_QUADRATIC_EDGE;
245 break;
246 case Rec4Face:
247 cellType = VTK_LINE;
248 break;
249 case Rec8Face:
250 cellType = VTK_QUADRATIC_EDGE;
251 break;
252 case Tet4Face:
253 cellType = VTK_TRIANGLE;
254 break;
255 case Tet10Face:
256 cellType = VTK_QUADRATIC_TRIANGLE;
257 break;
258 case Hex8Face:
259 cellType = VTK_QUAD;
260 break;
261 case Hex20Face:
262 cellType = VTK_QUADRATIC_QUAD;
263 break;
264 case Point1_Contact:
265 cellType = VTK_VERTEX;
266 break;
267 case Line2_Contact:
268 cellType = VTK_LINE;
269 break;
270 case Line3_Contact:
271 cellType = VTK_QUADRATIC_EDGE;
272 break;
273 case Tri3_Contact:
274 cellType = VTK_LINE;
275 break;
276 case Tri6_Contact:
277 cellType = VTK_QUADRATIC_TRIANGLE;
278 break;
279 case Rec4_Contact:
280 cellType = VTK_QUAD;
281 break;
282 case Rec8_Contact:
283 cellType = VTK_QUADRATIC_QUAD;
284 break;
285 case Line2Face_Contact:
286 cellType = VTK_VERTEX;
287 break;
288 case Line3Face_Contact:
289 cellType = VTK_VERTEX;
290 break;
291 case Tri3Face_Contact:
292 cellType = VTK_LINE;
293 break;
294 case Tri6Face_Contact:
295 cellType = VTK_QUADRATIC_EDGE;
296 break;
297 case Rec4Face_Contact:
298 cellType = VTK_LINE;
299 break;
300 case Rec8Face_Contact:
301 cellType = VTK_QUADRATIC_EDGE;
302 break;
303 case Tet4Face_Contact:
304 cellType = VTK_TRIANGLE;
305 break;
306 case Tet10Face_Contact:
307 cellType = VTK_QUADRATIC_TRIANGLE;
308 break;
309 case Hex8Face_Contact:
310 cellType = VTK_QUAD;
311 break;
312 case Hex20Face_Contact:
313 cellType = VTK_QUADRATIC_QUAD;
314 break;
315 default:
316 sprintf(error_msg, "saveVTK: Element type %s is not supported by VTK",elements->ReferenceElement->Type->Name);
317 Finley_setError(VALUE_ERROR,error_msg);
318 return;
319 }
320
321 switch(cellType) {
322 case VTK_VERTEX:
323 numVTKNodesPerElement = 1;
324 strcpy(elemTypeStr, "VTK_VERTEX");
325 break;
326 case VTK_LINE:
327 numVTKNodesPerElement = 2;
328 strcpy(elemTypeStr, "VTK_LINE");
329 break;
330 case VTK_TRIANGLE:
331 numVTKNodesPerElement = 3;
332 strcpy(elemTypeStr, "VTK_TRIANGLE");
333 break;
334 case VTK_QUAD:
335 numVTKNodesPerElement = 4;
336 strcpy(elemTypeStr, "VTK_QUAD");
337 break;
338 case VTK_TETRA:
339 numVTKNodesPerElement = 4;
340 strcpy(elemTypeStr, "VTK_TETRA");
341 break;
342 case VTK_HEXAHEDRON:
343 numVTKNodesPerElement = 8;
344 strcpy(elemTypeStr, "VTK_HEXAHEDRON");
345 break;
346 case VTK_QUADRATIC_EDGE:
347 numVTKNodesPerElement = 3;
348 strcpy(elemTypeStr, "VTK_QUADRATIC_EDGE");
349 break;
350 case VTK_QUADRATIC_TRIANGLE:
351 numVTKNodesPerElement = 6;
352 strcpy(elemTypeStr, "VTK_QUADRATIC_TRIANGLE");
353 break;
354 case VTK_QUADRATIC_QUAD:
355 numVTKNodesPerElement = 8;
356 strcpy(elemTypeStr, "VTK_QUADRATIC_QUAD");
357 break;
358 case VTK_QUADRATIC_TETRA:
359 numVTKNodesPerElement = 10;
360 strcpy(elemTypeStr, "VTK_QUADRATIC_TETRA");
361 break;
362 case VTK_QUADRATIC_HEXAHEDRON:
363 numVTKNodesPerElement = 20;
364 strcpy(elemTypeStr, "VTK_QUADRATIC_HEXAHEDRON");
365 break;
366 default:
367 sprintf(error_msg,"saveVTK: Cell type %d is not supported by VTK", cellType);
368 Finley_setError(VALUE_ERROR,error_msg);
369 return;
370 }
371
372 /* write out the DataArray element for the connectivity */
373 int NN = elements->ReferenceElement->Type->numNodes;
374 fprintf(fileHandle_p, "<Cells>\n");
375 fprintf(fileHandle_p, "<DataArray "
376 "Name=\"connectivity\" "
377 "type=\"Int32\" "
378 "format=\"ascii\">\n");
379 if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
380 for (i = 0; i < numCells; i++) {
381 for (j = 0; j < numVTKNodesPerElement; j++) {
382 j2=elements->ReferenceElement->Type->linearNodes[j];
383 fprintf(fileHandle_p,"%d ",mesh_p->Nodes->toReduced[elements->Nodes[INDEX2(j2, i, NN)]]);
384 }
385 fprintf(fileHandle_p, "\n");
386 }
387 } else if (VTK_QUADRATIC_HEXAHEDRON==cellType) {
388 for (i = 0; i < numCells; i++) {
389 fprintf(fileHandle_p,"%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d\n",
390 elements->Nodes[INDEX2(0, i, NN)],
391 elements->Nodes[INDEX2(1, i, NN)],
392 elements->Nodes[INDEX2(2, i, NN)],
393 elements->Nodes[INDEX2(3, i, NN)],
394 elements->Nodes[INDEX2(4, i, NN)],
395 elements->Nodes[INDEX2(5, i, NN)],
396 elements->Nodes[INDEX2(6, i, NN)],
397 elements->Nodes[INDEX2(7, i, NN)],
398 elements->Nodes[INDEX2(8, i, NN)],
399 elements->Nodes[INDEX2(9, i, NN)],
400 elements->Nodes[INDEX2(10, i, NN)],
401 elements->Nodes[INDEX2(11, i, NN)],
402 elements->Nodes[INDEX2(16, i, NN)],
403 elements->Nodes[INDEX2(17, i, NN)],
404 elements->Nodes[INDEX2(18, i, NN)],
405 elements->Nodes[INDEX2(19, i, NN)],
406 elements->Nodes[INDEX2(12, i, NN)],
407 elements->Nodes[INDEX2(13, i, NN)],
408 elements->Nodes[INDEX2(14, i, NN)],
409 elements->Nodes[INDEX2(15, i, NN)]);
410 }
411 } else if (numVTKNodesPerElement!=NN) {
412 for (i = 0; i < numCells; i++) {
413 for (j = 0; j < numVTKNodesPerElement; j++) {
414 j2=elements->ReferenceElement->Type->geoNodes[j];
415 fprintf(fileHandle_p,"%d ", elements->Nodes[INDEX2(j2, i, NN)]);
416 }
417 fprintf(fileHandle_p, "\n");
418 }
419 } else {
420 for (i = 0; i < numCells; i++) {
421 for (j = 0; j < numVTKNodesPerElement; j++) fprintf(fileHandle_p,"%d ", elements->Nodes[INDEX2(j, i, NN)]);
422 fprintf(fileHandle_p, "\n");
423 }
424 }
425 fprintf(fileHandle_p, "</DataArray>\n");
426
427 /* write out the DataArray element for the offsets */
428 fprintf(fileHandle_p, "<DataArray "
429 "Name=\"offsets\" "
430 "type=\"Int32\" "
431 "format=\"ascii\">\n");
432 for (i=numVTKNodesPerElement; i<=numCells*numVTKNodesPerElement; i+=numVTKNodesPerElement)
433 fprintf(fileHandle_p, "%d\n", i);
434 fprintf(fileHandle_p, "</DataArray>\n");
435
436 /* write out the DataArray element for the types */
437 fprintf(fileHandle_p, "<DataArray "
438 "Name=\"types\" "
439 "type=\"UInt8\" "
440 "format=\"ascii\">\n");
441 for (i=0; i<numCells; i++)
442 fprintf(fileHandle_p, "%d\n", cellType);
443 fprintf(fileHandle_p, "</DataArray>\n");
444
445 /* finish off the <Cells> element */
446 fprintf(fileHandle_p, "</Cells>\n");
447 }
448
449 /* data */
450 if (!isEmpty(data_p)) {
451 int rank = getDataPointRank(data_p);
452 int nComp = getDataPointSize(data_p);
453 /* barf if rank is greater than two */
454 if (rank > 2) {
455 sprintf(error_msg, "saveVTK: Vtk can't handle objects with rank greater than 2. object rank = %d\n", rank);
456 Finley_setError(VALUE_ERROR,error_msg);
457 return;
458 }
459 /* if the rank == 0: --> scalar data
460 * if the rank == 1: --> vector data
461 * if the rank == 2: --> tensor data
462 */
463 char dataNameStr[31], dataTypeStr[63];
464 int nCompReqd=1; /* the number of components required by vtk */
465 if (rank == 0) {
466 strcpy(dataNameStr, "scalar");
467 sprintf(dataTypeStr, "Scalars=\"%s\"", dataNameStr);
468 nCompReqd = 1;
469 }
470 else if (rank == 1) {
471 strcpy(dataNameStr, "vector");
472 sprintf(dataTypeStr, "Vectors=\"%s\"", dataNameStr);
473 nCompReqd = 3;
474 }
475 else if (rank == 2) {
476 strcpy(dataNameStr, "tensor");
477 sprintf(dataTypeStr, "Tensors=\"%s\"", dataNameStr);
478 nCompReqd = 9;
479 }
480 /* if have cell centred data then use CellData element,
481 * if have node centred data, then use PointData element
482 */
483 if (isCellCentered) {
484 /* now for the cell data */
485 fprintf(fileHandle_p, "<CellData %s>\n", dataTypeStr);
486 fprintf(fileHandle_p,
487 "<DataArray "
488 "Name=\"%s\" "
489 "type=\"Float32\" "
490 "NumberOfComponents=\"%d\" "
491 "format=\"ascii\">\n",
492 dataNameStr, nCompReqd);
493 int numPointsPerSample = elements->ReferenceElement->numQuadNodes;
494 if (numPointsPerSample) {
495 int shape = getDataPointShape(data_p, 0);
496 if (shape > 3) {
497 sprintf(error_msg,"saveVTK: shape should be 1, 2, or 3; I got %d\n", shape);
498 Finley_setError(VALUE_ERROR,error_msg);
499 return;
500 }
501 for (i=0; i<numCells; i++) {
502 values = getSampleData(data_p, i);
503 double sampleAvg[nComp];
504 for (k=0; k<nComp; k++) {
505 /* averaging over the number of points in the sample */
506 rtmp = 0.;
507 for (j=0; j<numPointsPerSample; j++)
508 rtmp += values[INDEX2(k,j,nComp)];
509 sampleAvg[k] = rtmp/numPointsPerSample;
510 }
511 /* if the number of required components is more than the number
512 * of actual components, pad with zeros
513 */
514 /* probably only need to get shape of first element */
515 /* write the data different ways for scalar, vector and tensor */
516 if (nCompReqd == 1) {
517 fprintf(fileHandle_p, " %e", sampleAvg[0]);
518 } else if (nCompReqd == 3) {
519 /* write out the data */
520 for (int m=0; m<shape; m++) {
521 fprintf(fileHandle_p, " %e", sampleAvg[m]);
522 }
523 /* pad with zeros */
524 for (int m=0; m<nCompReqd-shape; m++)
525 fprintf(fileHandle_p, " %e", 0.);
526 } else if (nCompReqd == 9) {
527 /* tensor data, so have a 3x3 matrix to output as a row
528 * of 9 data points */
529 int count = 0;
530 for (int m=0; m<shape; m++) {
531 for (int n=0; n<shape; n++) {
532 fprintf(fileHandle_p, " %e", sampleAvg[count]);
533 count++;
534 }
535 for (int n=0; n<3-shape; n++)
536 fprintf(fileHandle_p, " %e", 0.);
537 }
538 for (int m=0; m<3-shape; m++) {
539 for (int n=0; n<3; n++)
540 fprintf(fileHandle_p, " %e", 0.);
541 }
542 }
543 fprintf(fileHandle_p, "\n");
544 }
545 }
546 fprintf(fileHandle_p, "</DataArray>\n");
547 fprintf(fileHandle_p, "</CellData>\n");
548 } else {
549 /* now for the point data */
550 fprintf(fileHandle_p, "<PointData %s>\n", dataTypeStr);
551 fprintf(fileHandle_p, "<DataArray "
552 "Name=\"%s\" "
553 "type=\"Float32\" "
554 "NumberOfComponents=\"%d\" "
555 "format=\"ascii\">\n",
556 dataNameStr, nCompReqd);
557 /* write out the data */
558 /* if the number of required components is more than the number
559 * of actual components, pad with zeros
560 */
561 bool_t do_write=TRUE;
562 int shape = getDataPointShape(data_p, 0);
563 if (shape > 3) {
564 sprintf(error_msg,"shape should be 1, 2, or 3; I got %d\n", shape);
565 Finley_setError(VALUE_ERROR,error_msg);
566 return;
567 }
568 for (i=0; i<mesh_p->Nodes->numNodes; i++) {
569 switch (nodetype) {
570 case(FINLEY_DEGREES_OF_FREEDOM):
571 values = getSampleData(data_p,mesh_p->Nodes->degreeOfFreedom[i]);
572 break;
573 case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):
574 if (mesh_p->Nodes->toReduced[i]>=0) {
575 do_write=TRUE;
576 values = getSampleData(data_p,mesh_p->Nodes->reducedDegreeOfFreedom[i]);
577 } else {
578 do_write=FALSE;
579 }
580 break;
581 case(FINLEY_NODES):
582 values = getSampleData(data_p,i);
583 break;
584 }
585 /* write the data different ways for scalar, vector and tensor */
586 if (do_write) {
587 if (nCompReqd == 1) {
588 fprintf(fileHandle_p, " %e", values[0]);
589 }
590 else if (nCompReqd == 3) {
591 /* write out the data */
592 for (int m=0; m<shape; m++)
593 fprintf(fileHandle_p, " %e", values[m]);
594 /* pad with zeros */
595 for (int m=0; m<nCompReqd-shape; m++)
596 fprintf(fileHandle_p, " %e", 0.);
597 }
598 else if (nCompReqd == 9) {
599 /* tensor data, so have a 3x3 matrix to output as a row
600 * of 9 data points */
601 int count = 0;
602 for (int m=0; m<shape; m++) {
603 for (int n=0; n<shape; n++) {
604 fprintf(fileHandle_p, " %e", values[count]);
605 count++;
606 }
607 for (int n=0; n<3-shape; n++)
608 fprintf(fileHandle_p, " %e", 0.);
609 }
610 for (int m=0; m<3-shape; m++) {
611 for (int n=0; n<3; n++)
612 fprintf(fileHandle_p, " %e", 0.);
613 }
614 }
615 fprintf(fileHandle_p, "\n");
616 }
617 }
618 fprintf(fileHandle_p, "</DataArray>\n");
619 fprintf(fileHandle_p, "</PointData>\n");
620 }
621 }
622
623 /* finish off the piece */
624 fprintf(fileHandle_p, "</Piece>\n");
625
626 fprintf(fileHandle_p, "</UnstructuredGrid>\n");
627 /* write the xml footer */
628 fprintf(fileHandle_p, "</VTKFile>\n");
629 /* close the file */
630 fclose(fileHandle_p);
631 return;
632 }
633
634 /*
635 * Revision 1.6 2005/08/12 01:45:43 jgs
636 * erge of development branch dev-02 back to main trunk on 2005-08-12
637 *
638 * Revision 1.5.2.4 2005/09/09 08:15:17 gross
639 * bugs in vtk and dx writer fixed
640 *
641 * Revision 1.5.2.3 2005/09/08 08:28:39 gross
642 * some cleanup in savevtk
643 *
644 * Revision 1.5.2.2 2005/09/07 06:26:20 gross
645 * the solver from finley are put into the standalone package paso now
646 *
647 * Revision 1.5.2.1 2005/08/10 06:14:37 gross
648 * QUADRATIC HEXAHEDRON elements fixed
649 *
650 * Revision 1.5 2005/07/08 04:07:55 jgs
651 * Merge of development branch back to main trunk on 2005-07-08
652 *
653 * Revision 1.4 2005/05/06 04:26:15 jgs
654 * Merge of development branch back to main trunk on 2005-05-06
655 *
656 * Revision 1.1.2.7 2005/06/29 02:34:54 gross
657 * some changes towards 64 integers in finley
658 *
659 * Revision 1.1.2.6 2005/05/06 01:17:19 cochrane
660 * Fixed incorrect reporting of number of components in PointData arrays for
661 * vector data.
662 *
663 * Revision 1.1.2.5 2005/05/05 05:38:44 cochrane
664 * Improved formatting of VTK file output.
665 *
666 * Revision 1.1.2.4 2005/02/22 10:03:54 cochrane
667 * Implementation of writing of vtk xml files from finley. This function will
668 * require more testing, but on the cases that I have tried (and with the help
669 * of Lutz and mayavi), it looks like it's producing the correct output. Testing
670 * with more realistic data would be good. I'm at least confident enough to
671 * commit my changes.
672 *
673 * Revision 1.1.2.3 2005/02/17 05:53:26 gross
674 * some bug in saveDX fixed: in fact the bug was in
675 * DataC/getDataPointShape
676 *
677 * Revision 1.1.2.2 2005/02/10 01:34:22 cochrane
678 * Quick fix to make sure that saveVTK compiles so that finley is still buildable. Apologies to those this has affected.
679 *
680 * Revision 1.1.2.1 2005/02/09 06:53:15 cochrane
681 * Initial import to repository. This is the file to implement saving finley/escript meshes out to vtk formatted files. It is basically just a hack of the opendx equivalent, with a lot of the opendx stuff still in the file, so it doesn't actually work just yet, but it probably needs to be added to the cvs.
682 *
683 * Revision 1.1.1.1 2004/10/26 06:53:57 jgs
684 * initial import of project esys2
685 *
686 * Revision 1.1 2004/07/27 08:27:11 gross
687 * Finley: saveDX added: now it is possible to write data on boundary and contact elements
688 *
689 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26