/[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 687 - (show annotations)
Mon Mar 27 22:17:32 2006 UTC (13 years, 7 months ago) by gross
File MIME type: text/plain
File size: 27350 byte(s)
small fixes
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 /* Version: $Id$ */
22
23 /**************************************************************/
24
25 #include "Mesh.h"
26 #include "vtkCellType.h" /* copied from vtk source directory !!! */
27
28 /**************************************************************/
29
30 void Finley_Mesh_saveVTK(const char * filename_p, Finley_Mesh *mesh_p, const dim_t num_data,char* *names_p, escriptDataC* *data_pp) {
31 char error_msg[LenErrorMsg_MAX];
32 /* if there is no mesh we just return */
33 if (mesh_p==NULL) return;
34
35 int i, j, k, numVTKNodesPerElement,i_data,m, count, n, rank,shape, numPoints, cellType, numCells,
36 nDim, numPointsPerSample, nComp, nCompReqd;
37
38 index_t j2;
39 double* values, rtmp;
40 char elemTypeStr[32];
41
42 /* open the file and check handle */
43
44 FILE * fileHandle_p = fopen(filename_p, "w");
45 if (fileHandle_p==NULL) {
46 sprintf(error_msg, "saveVTK: File %s could not be opened for writing.", filename_p);
47 Finley_setError(IO_ERROR,error_msg);
48 return;
49 }
50 /* find the mesh type to be written */
51 int nodetype=FINLEY_DEGREES_OF_FREEDOM;
52 int elementtype=FINLEY_UNKNOWN;
53 bool_t isCellCentered[num_data],write_celldata=FALSE,write_pointdata=FALSE;
54 for (i_data=0;i_data<num_data;++i_data) {
55 if (! isEmpty(data_pp[i_data])) {
56 switch(getFunctionSpaceType(data_pp[i_data])) {
57 case FINLEY_DEGREES_OF_FREEDOM:
58 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
59 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
60 elementtype=FINLEY_ELEMENTS;
61 } else {
62 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
63 fclose(fileHandle_p);
64 return;
65 }
66 isCellCentered[i_data]=FALSE;
67 break;
68 case FINLEY_REDUCED_DEGREES_OF_FREEDOM:
69 nodetype = FINLEY_REDUCED_DEGREES_OF_FREEDOM;
70 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
71 elementtype=FINLEY_ELEMENTS;
72 } else {
73 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
74 fclose(fileHandle_p);
75 return;
76 }
77 isCellCentered[i_data]=FALSE;
78 break;
79 case FINLEY_NODES:
80 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
81 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
82 elementtype=FINLEY_ELEMENTS;
83 } else {
84 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
85 fclose(fileHandle_p);
86 return;
87 }
88 isCellCentered[i_data]=FALSE;
89 break;
90 case FINLEY_ELEMENTS:
91 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
92 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
93 elementtype=FINLEY_ELEMENTS;
94 } else {
95 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
96 fclose(fileHandle_p);
97 return;
98 }
99 isCellCentered[i_data]=TRUE;
100 break;
101 case FINLEY_FACE_ELEMENTS:
102 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
103 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_FACE_ELEMENTS) {
104 elementtype=FINLEY_FACE_ELEMENTS;
105 } else {
106 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
107 fclose(fileHandle_p);
108 return;
109 }
110 isCellCentered[i_data]=TRUE;
111 break;
112 case FINLEY_POINTS:
113 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
114 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_POINTS) {
115 elementtype=FINLEY_POINTS;
116 } else {
117 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
118 fclose(fileHandle_p);
119 return;
120 }
121 isCellCentered[i_data]=TRUE;
122 break;
123 case FINLEY_CONTACT_ELEMENTS_1:
124 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
125 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {
126 elementtype=FINLEY_CONTACT_ELEMENTS_1;
127 } else {
128 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
129 fclose(fileHandle_p);
130 return;
131 }
132 isCellCentered[i_data]=TRUE;
133 break;
134 case FINLEY_CONTACT_ELEMENTS_2:
135 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
136 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {
137 elementtype=FINLEY_CONTACT_ELEMENTS_1;
138 } else {
139 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
140 fclose(fileHandle_p);
141 return;
142 }
143 isCellCentered[i_data]=TRUE;
144 break;
145 default:
146 sprintf(error_msg,"saveVTK: Finley does not know anything about function space type %d",getFunctionSpaceType(data_pp[i_data]));
147 Finley_setError(TYPE_ERROR,error_msg);
148 fclose(fileHandle_p);
149 return;
150 }
151 if (isCellCentered[i_data]) {
152 write_celldata=TRUE;
153 } else {
154 write_pointdata=TRUE;
155 }
156 }
157 }
158 /* select nomber of points and the mesh component */
159 numPoints = mesh_p->Nodes->numNodes;
160 if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
161 numPoints = mesh_p->Nodes->reducedNumNodes;
162 } else {
163 numPoints = mesh_p->Nodes->numNodes;
164 }
165 if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;
166 Finley_ElementFile* elements=NULL;
167 switch(elementtype) {
168 case FINLEY_ELEMENTS:
169 elements=mesh_p->Elements;
170 break;
171 case FINLEY_FACE_ELEMENTS:
172 elements=mesh_p->FaceElements;
173 break;
174 case FINLEY_POINTS:
175 elements=mesh_p->Points;
176 break;
177 case FINLEY_CONTACT_ELEMENTS_1:
178 elements=mesh_p->ContactElements;
179 break;
180 }
181 if (elements==NULL) {
182 Finley_setError(SYSTEM_ERROR,"saveVTK: undefined element file");
183 fclose(fileHandle_p);
184 return;
185 }
186 /* map finley element type to VTK element type */
187 numCells = elements->numElements;
188 ElementTypeId TypeId;
189 if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
190 TypeId = elements->LinearReferenceElement->Type->TypeId;
191 } else {
192 TypeId = elements->ReferenceElement->Type->TypeId;
193 }
194
195 switch(TypeId) {
196 case Point1:
197 case Line2Face:
198 case Line3Face:
199 case Point1_Contact:
200 case Line2Face_Contact:
201 case Line3Face_Contact:
202 cellType = VTK_VERTEX;
203 numVTKNodesPerElement = 1;
204 strcpy(elemTypeStr, "VTK_VERTEX");
205 break;
206
207 case Line2:
208 case Tri3Face:
209 case Rec4Face:
210 case Line2_Contact:
211 case Tri3_Contact:
212 case Tri3Face_Contact:
213 case Rec4Face_Contact:
214 cellType = VTK_LINE;
215 numVTKNodesPerElement = 2;
216 strcpy(elemTypeStr, "VTK_LINE");
217 break;
218
219 case Tri3:
220 case Tet4Face:
221 case Tet4Face_Contact:
222 cellType = VTK_TRIANGLE;
223 numVTKNodesPerElement = 3;
224 strcpy(elemTypeStr, "VTK_TRIANGLE");
225 break;
226
227 case Rec4:
228 case Hex8Face:
229 case Rec4_Contact:
230 case Hex8Face_Contact:
231 cellType = VTK_QUAD;
232 numVTKNodesPerElement = 4;
233 strcpy(elemTypeStr, "VTK_QUAD");
234 break;
235
236 case Tet4:
237 cellType = VTK_TETRA;
238 numVTKNodesPerElement = 4;
239 strcpy(elemTypeStr, "VTK_TETRA");
240 break;
241
242 case Hex8:
243 cellType = VTK_HEXAHEDRON;
244 numVTKNodesPerElement = 8;
245 strcpy(elemTypeStr, "VTK_HEXAHEDRON");
246 break;
247
248 case Line3:
249 case Tri6Face:
250 case Rec8Face:
251 case Line3_Contact:
252 case Tri6Face_Contact:
253 case Rec8Face_Contact:
254 cellType = VTK_QUADRATIC_EDGE;
255 numVTKNodesPerElement = 3;
256 strcpy(elemTypeStr, "VTK_QUADRATIC_EDGE");
257 break;
258
259 case Tri6:
260 case Tet10Face:
261 case Tri6_Contact:
262 case Tet10Face_Contact:
263 cellType = VTK_QUADRATIC_TRIANGLE;
264 numVTKNodesPerElement = 6;
265 strcpy(elemTypeStr, "VTK_QUADRATIC_TRIANGLE");
266 break;
267
268 case Rec8:
269 case Hex20Face:
270 case Rec8_Contact:
271 case Hex20Face_Contact:
272 cellType = VTK_QUADRATIC_QUAD;
273 numVTKNodesPerElement = 8;
274 strcpy(elemTypeStr, "VTK_QUADRATIC_QUAD");
275 break;
276
277 case Tet10:
278 cellType = VTK_QUADRATIC_TETRA;
279 numVTKNodesPerElement = 10;
280 strcpy(elemTypeStr, "VTK_QUADRATIC_TETRA");
281 break;
282
283 case Hex20:
284 cellType = VTK_QUADRATIC_HEXAHEDRON;
285 numVTKNodesPerElement = 20;
286 strcpy(elemTypeStr, "VTK_QUADRATIC_HEXAHEDRON");
287 break;
288
289 default:
290 sprintf(error_msg, "saveVTK: Element type %s is not supported by VTK",elements->ReferenceElement->Type->Name);
291 Finley_setError(VALUE_ERROR,error_msg);
292 fclose(fileHandle_p);
293 return;
294 }
295 /* xml header */
296 fprintf(fileHandle_p, "<?xml version=\"1.0\"?>\n");
297 fprintf(fileHandle_p, "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\">\n");
298
299 /* finley uses an unstructured mesh, so UnstructuredGrid *should* work */
300 fprintf(fileHandle_p, "<UnstructuredGrid>\n");
301
302 /* is there only one "piece" to the data?? */
303 fprintf(fileHandle_p, "<Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n",numPoints, numCells);
304 /* now for the points; equivalent to positions section in saveDX() */
305 /* "The points element explicitly defines coordinates for each point
306 * individually. It contains one DataArray element describing an array
307 * with three components per value, each specifying the coordinates of one
308 * point" - from Vtk User's Guide
309 */
310 fprintf(fileHandle_p, "<Points>\n");
311 /*
312 * the reason for this if statement is explained in the long comment below
313 */
314 nDim = mesh_p->Nodes->numDim;
315 fprintf(fileHandle_p, "<DataArray NumberOfComponents=\"%d\" type=\"Float32\" format=\"ascii\">\n",MAX(3,nDim));
316 /* vtk/mayavi doesn't like 2D data, it likes 3D data with a degenerate
317 * third dimension to handle 2D data (like a sheet of paper). So, if
318 * nDim is 2, we have to append zeros to the array to get this third
319 * dimension, and keep the visualisers happy.
320 * Indeed, if nDim is less than 3, must pad all empty dimensions, so
321 * that the total number of dims is 3.
322 */
323 if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
324 for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
325 if (mesh_p->Nodes->toReduced[i]>=0) {
326 for (j = 0; j < nDim; j++) fprintf(fileHandle_p, " %e",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);
327 for (k=0; k<3-nDim; k++) fprintf(fileHandle_p, " %e",0.);
328 fprintf(fileHandle_p, "\n");
329 }
330 }
331 } else {
332 for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
333 for (j = 0; j < nDim; j++) fprintf(fileHandle_p, " %e",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);
334 for (k=0; k<3-nDim; k++) fprintf(fileHandle_p, " %e",0.);
335 fprintf(fileHandle_p, "\n");
336 }
337 }
338 fprintf(fileHandle_p, "</DataArray>\n");
339 fprintf(fileHandle_p, "</Points>\n");
340
341 /* write out the DataArray element for the connectivity */
342
343 int NN = elements->ReferenceElement->Type->numNodes;
344 fprintf(fileHandle_p, "<Cells>\n");
345 fprintf(fileHandle_p, "<DataArray Name=\"connectivity\" type=\"Int32\" format=\"ascii\">\n");
346
347 if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
348 for (i = 0; i < numCells; i++) {
349 for (j = 0; j < numVTKNodesPerElement; j++)
350 fprintf(fileHandle_p,"%d ",mesh_p->Nodes->toReduced[elements->Nodes[INDEX2(elements->ReferenceElement->Type->linearNodes[j], i, NN)]]);
351 fprintf(fileHandle_p, "\n");
352 }
353 } else if (VTK_QUADRATIC_HEXAHEDRON==cellType) {
354 for (i = 0; i < numCells; i++) {
355 fprintf(fileHandle_p,"%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d\n",
356 elements->Nodes[INDEX2(0, i, NN)],
357 elements->Nodes[INDEX2(1, i, NN)],
358 elements->Nodes[INDEX2(2, i, NN)],
359 elements->Nodes[INDEX2(3, i, NN)],
360 elements->Nodes[INDEX2(4, i, NN)],
361 elements->Nodes[INDEX2(5, i, NN)],
362 elements->Nodes[INDEX2(6, i, NN)],
363 elements->Nodes[INDEX2(7, i, NN)],
364 elements->Nodes[INDEX2(8, i, NN)],
365 elements->Nodes[INDEX2(9, i, NN)],
366 elements->Nodes[INDEX2(10, i, NN)],
367 elements->Nodes[INDEX2(11, i, NN)],
368 elements->Nodes[INDEX2(16, i, NN)],
369 elements->Nodes[INDEX2(17, i, NN)],
370 elements->Nodes[INDEX2(18, i, NN)],
371 elements->Nodes[INDEX2(19, i, NN)],
372 elements->Nodes[INDEX2(12, i, NN)],
373 elements->Nodes[INDEX2(13, i, NN)],
374 elements->Nodes[INDEX2(14, i, NN)],
375 elements->Nodes[INDEX2(15, i, NN)]);
376 }
377 } else if (numVTKNodesPerElement!=NN) {
378 for (i = 0; i < numCells; i++) {
379 for (j = 0; j < numVTKNodesPerElement; j++) fprintf(fileHandle_p,"%d ", elements->Nodes[INDEX2(elements->ReferenceElement->Type->geoNodes[j], i, NN)]);
380 fprintf(fileHandle_p, "\n");
381 }
382 } else {
383 for (i = 0; i < numCells; i++) {
384 for (j = 0; j < numVTKNodesPerElement; j++) fprintf(fileHandle_p,"%d ", elements->Nodes[INDEX2(j, i, NN)]);
385 fprintf(fileHandle_p, "\n");
386 }
387 }
388 fprintf(fileHandle_p, "</DataArray>\n");
389
390 /* write out the DataArray element for the offsets */
391 fprintf(fileHandle_p, "<DataArray Name=\"offsets\" type=\"Int32\" format=\"ascii\">\n");
392 for (i=numVTKNodesPerElement; i<=numCells*numVTKNodesPerElement; i+=numVTKNodesPerElement) fprintf(fileHandle_p, "%d\n", i);
393 fprintf(fileHandle_p, "</DataArray>\n");
394
395 /* write out the DataArray element for the types */
396 fprintf(fileHandle_p, "<DataArray Name=\"types\" type=\"UInt8\" format=\"ascii\">\n");
397 for (i=0; i<numCells; i++) fprintf(fileHandle_p, "%d\n", cellType);
398 fprintf(fileHandle_p, "</DataArray>\n");
399
400 /* finish off the <Cells> element */
401 fprintf(fileHandle_p, "</Cells>\n");
402
403 /* cell data */
404 if (write_celldata) {
405 /* mark the active data arrays */
406 bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
407 fprintf(fileHandle_p, "<CellData");
408 for (i_data =0 ;i_data<num_data;++i_data) {
409 if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data]) {
410 /* if the rank == 0: --> scalar data
411 * if the rank == 1: --> vector data
412 * if the rank == 2: --> tensor data
413 */
414 switch(getDataPointRank(data_pp[i_data])) {
415 case 0:
416 if (! set_scalar) {
417 fprintf(fileHandle_p," Scalars=\"%s\"",names_p[i_data]);
418 set_scalar=TRUE;
419 }
420 break;
421 case 1:
422 if (! set_vector) {
423 fprintf(fileHandle_p," Vectors=\"%s\"",names_p[i_data]);
424 set_vector=TRUE;
425 }
426 break;
427 case 2:
428 if (! set_tensor) {
429 fprintf(fileHandle_p," Tensors=\"%s\"",names_p[i_data]);
430 set_tensor=TRUE;
431 }
432 break;
433 default:
434 sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);
435 Finley_setError(VALUE_ERROR,error_msg);
436 fclose(fileHandle_p);
437 return;
438 }
439 }
440 }
441 fprintf(fileHandle_p, ">\n");
442 /* write the arrays */
443 for (i_data =0 ;i_data<num_data;++i_data) {
444 if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data]) {
445 numPointsPerSample = elements->ReferenceElement->numQuadNodes;
446 rank = getDataPointRank(data_pp[i_data]);
447 nComp = getDataPointSize(data_pp[i_data]);
448 nCompReqd=1; /* the number of components required by vtk */
449 shape=0;
450 if (rank == 0) {
451 nCompReqd = 1;
452 } else if (rank == 1) {
453 shape=getDataPointShape(data_pp[i_data], 0);
454 if (shape>3) {
455 Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");
456 fclose(fileHandle_p);
457 return;
458 }
459 nCompReqd = 3;
460 } else {
461 shape=getDataPointShape(data_pp[i_data], 0);
462 if (shape>3 || shape != getDataPointShape(data_pp[i_data], 1)) {
463 Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
464 fclose(fileHandle_p);
465 return;
466 }
467 nCompReqd = 9;
468 }
469 fprintf(fileHandle_p, "<DataArray Name=\"%s\" type=\"Float32\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);
470
471 double sampleAvg[nComp];
472 for (i=0; i<numCells; i++) {
473 values = getSampleData(data_pp[i_data], i);
474 /* averaging over the number of points in the sample */
475 for (k=0; k<nComp; k++) {
476 rtmp = 0.;
477 for (j=0; j<numPointsPerSample; j++) rtmp += values[INDEX2(k,j,nComp)];
478 sampleAvg[k] = rtmp/numPointsPerSample;
479 }
480 /* if the number of required components is more than the number
481 * of actual components, pad with zeros
482 */
483 /* probably only need to get shape of first element */
484 /* write the data different ways for scalar, vector and tensor */
485 if (nCompReqd == 1) {
486 fprintf(fileHandle_p, " %e", sampleAvg[0]);
487 } else if (nCompReqd == 3) {
488 /* write out the data */
489 for (m=0; m<shape; m++) fprintf(fileHandle_p, " %e", sampleAvg[m]);
490 for (m=0; m<nCompReqd-shape; m++) fprintf(fileHandle_p, " %e", 0.);
491 } else if (nCompReqd == 9) {
492 /* tensor data, so have a 3x3 matrix to output as a row
493 * of 9 data points */
494 count = 0;
495 for (m=0; m<shape; m++) {
496 for (n=0; n<shape; n++) {
497 fprintf(fileHandle_p, " %e", sampleAvg[count]);
498 count++;
499 }
500 for (n=0; n<3-shape; n++) fprintf(fileHandle_p, " %e", 0.);
501 }
502 for (m=0; m<3-shape; m++)
503 for (n=0; n<3; n++) fprintf(fileHandle_p, " %e", 0.);
504 }
505 fprintf(fileHandle_p, "\n");
506 }
507 fprintf(fileHandle_p, "</DataArray>\n");
508 }
509 }
510 fprintf(fileHandle_p, "</CellData>\n");
511 }
512 /* point data */
513 if (write_pointdata) {
514 /* mark the active data arrays */
515 bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
516 fprintf(fileHandle_p, "<PointData");
517 for (i_data =0 ;i_data<num_data;++i_data) {
518 if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data]) {
519 /* if the rank == 0: --> scalar data
520 * if the rank == 1: --> vector data
521 * if the rank == 2: --> tensor data
522 */
523 switch(getDataPointRank(data_pp[i_data])) {
524 case 0:
525 if (! set_scalar) {
526 fprintf(fileHandle_p," Scalars=\"%s\"",names_p[i_data]);
527 set_scalar=TRUE;
528 }
529 break;
530 case 1:
531 if (! set_vector) {
532 fprintf(fileHandle_p," Vectors=\"%s\"",names_p[i_data]);
533 set_vector=TRUE;
534 }
535 break;
536 case 2:
537 if (! set_tensor) {
538 fprintf(fileHandle_p," Tensors=\"%s\"",names_p[i_data]);
539 set_tensor=TRUE;
540 }
541 break;
542 default:
543 sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);
544 Finley_setError(VALUE_ERROR,error_msg);
545 fclose(fileHandle_p);
546 return;
547 }
548 }
549 }
550 fprintf(fileHandle_p, ">\n");
551 /* write the arrays */
552 for (i_data =0 ;i_data<num_data;++i_data) {
553 if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data]) {
554 numPointsPerSample = elements->ReferenceElement->numQuadNodes;
555 rank = getDataPointRank(data_pp[i_data]);
556 nComp = getDataPointSize(data_pp[i_data]);
557 nCompReqd=1; /* the number of components required by vtk */
558 shape=0;
559 if (rank == 0) {
560 nCompReqd = 1;
561 } else if (rank == 1) {
562 shape=getDataPointShape(data_pp[i_data], 0);
563 if (shape>3) {
564 Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");
565 fclose(fileHandle_p);
566 return;
567 }
568 nCompReqd = 3;
569 } else {
570 shape=getDataPointShape(data_pp[i_data], 0);
571 if (shape>3 || shape != getDataPointShape(data_pp[i_data], 1)) {
572 Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
573 fclose(fileHandle_p);
574 return;
575 }
576 nCompReqd = 9;
577 }
578 fprintf(fileHandle_p, "<DataArray Name=\"%s\" type=\"Float32\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);
579 /* write out the data */
580 /* if the number of required components is more than the number
581 * of actual components, pad with zeros
582 */
583 bool_t do_write=TRUE;
584 for (i=0; i<mesh_p->Nodes->numNodes; i++) {
585 if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
586 if (mesh_p->Nodes->toReduced[i]>=0) {
587 switch(getFunctionSpaceType(data_pp[i_data])) {
588 case FINLEY_DEGREES_OF_FREEDOM:
589 values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);
590 break;
591 case FINLEY_REDUCED_DEGREES_OF_FREEDOM:
592 values = getSampleData(data_pp[i_data],mesh_p->Nodes->reducedDegreeOfFreedom[i]);
593 break;
594 case FINLEY_NODES:
595 values = getSampleData(data_pp[i_data],i);
596 break;
597 }
598 do_write=TRUE;
599 } else {
600 do_write=FALSE;
601 }
602 } else {
603 do_write=TRUE;
604 switch(getFunctionSpaceType(data_pp[i_data])) {
605 case FINLEY_DEGREES_OF_FREEDOM:
606 values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);
607 break;
608 case FINLEY_NODES:
609 values = getSampleData(data_pp[i_data],i);
610 break;
611 }
612 }
613 /* write the data different ways for scalar, vector and tensor */
614 if (do_write) {
615 if (nCompReqd == 1) {
616 fprintf(fileHandle_p, " %e", values[0]);
617 } else if (nCompReqd == 3) {
618 for (m=0; m<shape; m++) fprintf(fileHandle_p, " %e", values[m]);
619 for (m=0; m<nCompReqd-shape; m++) fprintf(fileHandle_p, " %e", 0.);
620 } else if (nCompReqd == 9) {
621 /* tensor data, so have a 3x3 matrix to output as a row
622 * of 9 data points */
623 count = 0;
624 for (m=0; m<shape; m++) {
625 for (n=0; n<shape; n++) {
626 fprintf(fileHandle_p, " %e", values[count]);
627 count++;
628 }
629 for (n=0; n<3-shape; n++) fprintf(fileHandle_p, " %e", 0.);
630 }
631 for (m=0; m<3-shape; m++)
632 for (n=0; n<3; n++) fprintf(fileHandle_p, " %e", 0.);
633 }
634 fprintf(fileHandle_p, "\n");
635 }
636 }
637 fprintf(fileHandle_p, "</DataArray>\n");
638 }
639 }
640 fprintf(fileHandle_p, "</PointData>\n");
641 }
642 /* finish off the piece */
643 fprintf(fileHandle_p, "</Piece>\n");
644
645 fprintf(fileHandle_p, "</UnstructuredGrid>\n");
646 /* write the xml footer */
647 fprintf(fileHandle_p, "</VTKFile>\n");
648 /* close the file */
649 fclose(fileHandle_p);
650 return;
651 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26