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

Contents of /branches/RW_WIN32/finley/src/finleyC/Mesh_saveVTK.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 210 - (show annotations)
Wed Nov 23 09:54:02 2005 UTC (14 years, 11 months ago) by robwdcock
File MIME type: text/plain
File size: 28980 byte(s)
PARTIAL WIN32 BUILD SYSTEM AND PORT
+ All libraries now build under new build system. No unit test routines yet
+ Reversed some incorrect refactorings in Bruce.cpp

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26