/[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 616 - (show annotations)
Wed Mar 22 02:46:56 2006 UTC (13 years, 3 months ago) by elspeth
File MIME type: text/plain
File size: 27376 byte(s)
Copyright added to more source files.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26