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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 272 - (hide annotations)
Wed Nov 30 03:30:30 2005 UTC (16 years, 8 months ago) by gross
File MIME type: text/plain
File size: 27830 byte(s) 
cast to float added. this avoids problems with very small coodinates.
1 jgs 150 /*
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 jgs 110
15 jgs 150
16 jgs 110 /**************************************************************/
17    
18     /* writes data and mesh in a vtk file */
19    
20     /**************************************************************/
21    
22     /* Author: Paul Cochrane, cochrane@esscc.uq.edu.au */
23 jgs 150 /* Version: $Id$ */
24 jgs 110
25     /**************************************************************/
26    
27     #include "Mesh.h"
28 jgs 113 #include "vtkCellType.h" /* copied from vtk source directory !!! */
29 jgs 110
30 jgs 150 /**************************************************************/
31    
32 jgs 153 void Finley_Mesh_saveVTK(const char * filename_p, Finley_Mesh *mesh_p, const dim_t num_data,char* *names_p, escriptDataC* *data_pp) {
33 jgs 150 char error_msg[LenErrorMsg_MAX];
34 jgs 110 /* if there is no mesh we just return */
35     if (mesh_p==NULL) return;
36 jgs 153
37     int i, j, k, numVTKNodesPerElement,i_data;
38 jgs 150 index_t j2;
39 jgs 113 double* values, rtmp;
40 jgs 150
41 jgs 153 /* 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 jgs 150 }
49 jgs 153 /* 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 jgs 113 } else {
162 jgs 150 numPoints = mesh_p->Nodes->numNodes;
163 jgs 153 }
164     if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;
165     Finley_ElementFile* elements=NULL;
166     switch(elementtype) {
167     case FINLEY_ELEMENTS:
168 jgs 113 elements=mesh_p->Elements;
169     break;
170 jgs 153 case FINLEY_FACE_ELEMENTS:
171 jgs 113 elements=mesh_p->FaceElements;
172     break;
173 jgs 153 case FINLEY_POINTS:
174 jgs 113 elements=mesh_p->Points;
175     break;
176 jgs 153 case FINLEY_CONTACT_ELEMENTS_1:
177 jgs 113 elements=mesh_p->ContactElements;
178     break;
179     }
180 jgs 153 if (elements==NULL) {
181     Finley_setError(SYSTEM_ERROR,"saveVTK: undefined element file");
182     fclose(fileHandle_p);
183     return;
184 jgs 113 }
185 jgs 153 /* map finley element type to VTK element type */
186 jgs 113 int numCells = elements->numElements;
187 jgs 153 int cellType;
188     ElementTypeId TypeId;
189     char elemTypeStr[32];
190     if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
191     TypeId = elements->LinearReferenceElement->Type->TypeId;
192 jgs 113 } else {
193 jgs 153 TypeId = elements->ReferenceElement->Type->TypeId;
194 jgs 113 }
195 jgs 110
196 jgs 153 switch(TypeId) {
197 jgs 113 case Point1:
198     case Line2Face:
199     case Line3Face:
200 jgs 153 case Point1_Contact:
201     case Line2Face_Contact:
202     case Line3Face_Contact:
203 jgs 113 cellType = VTK_VERTEX;
204 jgs 153 numVTKNodesPerElement = 1;
205     strcpy(elemTypeStr, "VTK_VERTEX");
206 jgs 113 break;
207 jgs 153
208     case Line2:
209 jgs 113 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 jgs 153 numVTKNodesPerElement = 2;
217     strcpy(elemTypeStr, "VTK_LINE");
218 jgs 113 break;
219 jgs 153
220     case Tri3:
221     case Tet4Face:
222 jgs 113 case Tet4Face_Contact:
223     cellType = VTK_TRIANGLE;
224 jgs 153 numVTKNodesPerElement = 3;
225     strcpy(elemTypeStr, "VTK_TRIANGLE");
226 jgs 113 break;
227 jgs 153
228     case Rec4:
229     case Hex8Face:
230     case Rec4_Contact:
231 jgs 113 case Hex8Face_Contact:
232     cellType = VTK_QUAD;
233     numVTKNodesPerElement = 4;
234     strcpy(elemTypeStr, "VTK_QUAD");
235     break;
236 jgs 153
237     case Tet4:
238     cellType = VTK_TETRA;
239 jgs 113 numVTKNodesPerElement = 4;
240     strcpy(elemTypeStr, "VTK_TETRA");
241     break;
242 jgs 153
243     case Hex8:
244     cellType = VTK_HEXAHEDRON;
245 jgs 113 numVTKNodesPerElement = 8;
246     strcpy(elemTypeStr, "VTK_HEXAHEDRON");
247     break;
248 jgs 153
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 jgs 113 numVTKNodesPerElement = 3;
257     strcpy(elemTypeStr, "VTK_QUADRATIC_EDGE");
258     break;
259 jgs 153
260     case Tri6:
261     case Tet10Face:
262     case Tri6_Contact:
263     case Tet10Face_Contact:
264     cellType = VTK_QUADRATIC_TRIANGLE;
265 jgs 113 numVTKNodesPerElement = 6;
266     strcpy(elemTypeStr, "VTK_QUADRATIC_TRIANGLE");
267     break;
268 jgs 153
269     case Rec8:
270     case Hex20Face:
271     case Rec8_Contact:
272     case Hex20Face_Contact:
273     cellType = VTK_QUADRATIC_QUAD;
274 jgs 113 numVTKNodesPerElement = 8;
275     strcpy(elemTypeStr, "VTK_QUADRATIC_QUAD");
276     break;
277 jgs 153
278     case Tet10:
279     cellType = VTK_QUADRATIC_TETRA;
280 jgs 113 numVTKNodesPerElement = 10;
281     strcpy(elemTypeStr, "VTK_QUADRATIC_TETRA");
282     break;
283 jgs 153
284     case Hex20:
285     cellType = VTK_QUADRATIC_HEXAHEDRON;
286 jgs 113 numVTKNodesPerElement = 20;
287     strcpy(elemTypeStr, "VTK_QUADRATIC_HEXAHEDRON");
288     break;
289 jgs 153
290     default:
291     sprintf(error_msg, "saveVTK: Element type %s is not supported by VTK",elements->ReferenceElement->Type->Name);
292 jgs 150 Finley_setError(VALUE_ERROR,error_msg);
293 jgs 153 fclose(fileHandle_p);
294 jgs 113 return;
295 jgs 153 }
296     /* xml header */
297     fprintf(fileHandle_p, "<?xml version=\"1.0\"?>\n");
298     fprintf(fileHandle_p, "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\">\n");
299 jgs 110
300 jgs 153 /* 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 gross 272 for (j = 0; j < nDim; j++) fprintf(fileHandle_p, " %e",(float) (mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]));
328     for (k=0; k<3-nDim; k++) fprintf(fileHandle_p, " %e",(float) 0.);
329 jgs 150 fprintf(fileHandle_p, "\n");
330 jgs 147 }
331 jgs 150 }
332 jgs 153 } else {
333     for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
334 gross 272 for (j = 0; j < nDim; j++) fprintf(fileHandle_p, " %e",(float) (mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]));
335     for (k=0; k<3-nDim; k++) fprintf(fileHandle_p, " %e",(float) 0.);
336 jgs 153 fprintf(fileHandle_p, "\n");
337     }
338     }
339     fprintf(fileHandle_p, "</DataArray>\n");
340     fprintf(fileHandle_p, "</Points>\n");
341 jgs 113
342 jgs 153 /* write out the DataArray element for the connectivity */
343 jgs 113
344 jgs 153 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 jgs 113
348 jgs 153 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 jgs 110
391 jgs 153 /* 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 jgs 150 }
441 jgs 153 }
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 gross 271 fprintf(fileHandle_p, " %e", (float) sampleAvg[0]);
488 jgs 153 } else if (nCompReqd == 3) {
489     /* write out the data */
490 gross 271 for (int m=0; m<shape; m++) fprintf(fileHandle_p, " %e", (float) sampleAvg[m]);
491     for (int m=0; m<nCompReqd-shape; m++) fprintf(fileHandle_p, " %e", (float) 0.);
492 jgs 153 } 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 gross 271 fprintf(fileHandle_p, " %e", (float) sampleAvg[count]);
499 jgs 153 count++;
500     }
501 gross 271 for (int n=0; n<3-shape; n++) fprintf(fileHandle_p, " %e", (float) 0.);
502 jgs 153 }
503     for (int m=0; m<3-shape; m++)
504 gross 271 for (int n=0; n<3; n++) fprintf(fileHandle_p, " %e", (float) 0.);
505 jgs 153 }
506     fprintf(fileHandle_p, "\n");
507     }
508     fprintf(fileHandle_p, "</DataArray>\n");
509 jgs 150 }
510 jgs 153 }
511     fprintf(fileHandle_p, "</CellData>\n");
512 jgs 110 }
513 jgs 153 /* 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 gross 271 fprintf(fileHandle_p, " %e", (float) values[0]);
618 jgs 153 } else if (nCompReqd == 3) {
619 gross 271 for (int m=0; m<shape; m++) fprintf(fileHandle_p, " %e", (float) values[m]);
620     for (int m=0; m<nCompReqd-shape; m++) fprintf(fileHandle_p, " %e", (float) 0.);
621 jgs 153 } 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 gross 271 fprintf(fileHandle_p, " %e", (float) values[count]);
628 jgs 153 count++;
629     }
630 gross 271 for (int n=0; n<3-shape; n++) fprintf(fileHandle_p, " %e", (float) 0.);
631 jgs 153 }
632     for (int m=0; m<3-shape; m++)
633 gross 271 for (int n=0; n<3; n++) fprintf(fileHandle_p, " %e", (float) 0.);
634 jgs 153 }
635     fprintf(fileHandle_p, "\n");
636     }
637     }
638     fprintf(fileHandle_p, "</DataArray>\n");
639     }
640     }
641     fprintf(fileHandle_p, "</PointData>\n");
642     }
643 jgs 113 /* finish off the piece */
644     fprintf(fileHandle_p, "</Piece>\n");
645 jgs 110
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