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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26