/[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 471 - (show annotations)
Fri Jan 27 01:33:02 2006 UTC (13 years, 3 months ago) by jgs
File MIME type: text/plain
File size: 27698 byte(s)
reorganise finley src tree to remove inc dir and src/finley directory

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26