/[escript]/branches/domexper/dudley/src/Mesh_saveVTK.c
ViewVC logotype

Annotation of /branches/domexper/dudley/src/Mesh_saveVTK.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 616 - (hide annotations)
Wed Mar 22 02:46:56 2006 UTC (13 years, 2 months ago) by elspeth
Original Path: trunk/finley/src/Mesh_saveVTK.c
File MIME type: text/plain
File size: 27376 byte(s)
Copyright added to more source files.

1 jgs 150 /*
2 elspeth 616 ************************************************************
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 jgs 150 */
12 jgs 110
13 jgs 150
14 jgs 110 /**************************************************************/
15    
16     /* writes data and mesh in a vtk file */
17    
18     /**************************************************************/
19    
20     /* Author: Paul Cochrane, cochrane@esscc.uq.edu.au */
21 jgs 150 /* Version: $Id$ */
22 jgs 110
23     /**************************************************************/
24    
25     #include "Mesh.h"
26 jgs 113 #include "vtkCellType.h" /* copied from vtk source directory !!! */
27 jgs 110
28 jgs 150 /**************************************************************/
29    
30 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) {
31 jgs 150 char error_msg[LenErrorMsg_MAX];
32 jgs 110 /* if there is no mesh we just return */
33     if (mesh_p==NULL) return;
34 jgs 153
35     int i, j, k, numVTKNodesPerElement,i_data;
36 jgs 150 index_t j2;
37 jgs 113 double* values, rtmp;
38 jgs 150
39 jgs 153 /* 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 jgs 150 }
47 jgs 153 /* 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 jgs 113 } else {
160 jgs 150 numPoints = mesh_p->Nodes->numNodes;
161 jgs 153 }
162     if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;
163     Finley_ElementFile* elements=NULL;
164     switch(elementtype) {
165     case FINLEY_ELEMENTS:
166 jgs 113 elements=mesh_p->Elements;
167     break;
168 jgs 153 case FINLEY_FACE_ELEMENTS:
169 jgs 113 elements=mesh_p->FaceElements;
170     break;
171 jgs 153 case FINLEY_POINTS:
172 jgs 113 elements=mesh_p->Points;
173     break;
174 jgs 153 case FINLEY_CONTACT_ELEMENTS_1:
175 jgs 113 elements=mesh_p->ContactElements;
176     break;
177     }
178 jgs 153 if (elements==NULL) {
179     Finley_setError(SYSTEM_ERROR,"saveVTK: undefined element file");
180     fclose(fileHandle_p);
181     return;
182 jgs 113 }
183 jgs 153 /* map finley element type to VTK element type */
184 jgs 113 int numCells = elements->numElements;
185 jgs 153 int cellType;
186     ElementTypeId TypeId;
187     char elemTypeStr[32];
188     if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
189     TypeId = elements->LinearReferenceElement->Type->TypeId;
190 jgs 113 } else {
191 jgs 153 TypeId = elements->ReferenceElement->Type->TypeId;
192 jgs 113 }
193 jgs 110
194 jgs 153 switch(TypeId) {
195 jgs 113 case Point1:
196     case Line2Face:
197     case Line3Face:
198 jgs 153 case Point1_Contact:
199     case Line2Face_Contact:
200     case Line3Face_Contact:
201 jgs 113 cellType = VTK_VERTEX;
202 jgs 153 numVTKNodesPerElement = 1;
203     strcpy(elemTypeStr, "VTK_VERTEX");
204 jgs 113 break;
205 jgs 153
206     case Line2:
207 jgs 113 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 jgs 153 numVTKNodesPerElement = 2;
215     strcpy(elemTypeStr, "VTK_LINE");
216 jgs 113 break;
217 jgs 153
218     case Tri3:
219     case Tet4Face:
220 jgs 113 case Tet4Face_Contact:
221     cellType = VTK_TRIANGLE;
222 jgs 153 numVTKNodesPerElement = 3;
223     strcpy(elemTypeStr, "VTK_TRIANGLE");
224 jgs 113 break;
225 jgs 153
226     case Rec4:
227     case Hex8Face:
228     case Rec4_Contact:
229 jgs 113 case Hex8Face_Contact:
230     cellType = VTK_QUAD;
231     numVTKNodesPerElement = 4;
232     strcpy(elemTypeStr, "VTK_QUAD");
233     break;
234 jgs 153
235     case Tet4:
236     cellType = VTK_TETRA;
237 jgs 113 numVTKNodesPerElement = 4;
238     strcpy(elemTypeStr, "VTK_TETRA");
239     break;
240 jgs 153
241     case Hex8:
242     cellType = VTK_HEXAHEDRON;
243 jgs 113 numVTKNodesPerElement = 8;
244     strcpy(elemTypeStr, "VTK_HEXAHEDRON");
245     break;
246 jgs 153
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 jgs 113 numVTKNodesPerElement = 3;
255     strcpy(elemTypeStr, "VTK_QUADRATIC_EDGE");
256     break;
257 jgs 153
258     case Tri6:
259     case Tet10Face:
260     case Tri6_Contact:
261     case Tet10Face_Contact:
262     cellType = VTK_QUADRATIC_TRIANGLE;
263 jgs 113 numVTKNodesPerElement = 6;
264     strcpy(elemTypeStr, "VTK_QUADRATIC_TRIANGLE");
265     break;
266 jgs 153
267     case Rec8:
268     case Hex20Face:
269     case Rec8_Contact:
270     case Hex20Face_Contact:
271     cellType = VTK_QUADRATIC_QUAD;
272 jgs 113 numVTKNodesPerElement = 8;
273     strcpy(elemTypeStr, "VTK_QUADRATIC_QUAD");
274     break;
275 jgs 153
276     case Tet10:
277     cellType = VTK_QUADRATIC_TETRA;
278 jgs 113 numVTKNodesPerElement = 10;
279     strcpy(elemTypeStr, "VTK_QUADRATIC_TETRA");
280     break;
281 jgs 153
282     case Hex20:
283     cellType = VTK_QUADRATIC_HEXAHEDRON;
284 jgs 113 numVTKNodesPerElement = 20;
285     strcpy(elemTypeStr, "VTK_QUADRATIC_HEXAHEDRON");
286     break;
287 jgs 153
288     default:
289     sprintf(error_msg, "saveVTK: Element type %s is not supported by VTK",elements->ReferenceElement->Type->Name);
290 jgs 150 Finley_setError(VALUE_ERROR,error_msg);
291 jgs 153 fclose(fileHandle_p);
292 jgs 113 return;
293 jgs 153 }
294     /* xml header */
295     fprintf(fileHandle_p, "<?xml version=\"1.0\"?>\n");
296     fprintf(fileHandle_p, "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\">\n");
297 jgs 110
298 jgs 153 /* 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 jgs 150 fprintf(fileHandle_p, "\n");
328 jgs 147 }
329 jgs 150 }
330 jgs 153 } 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 jgs 113
340 jgs 153 /* write out the DataArray element for the connectivity */
341 jgs 113
342 jgs 153 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 jgs 113
346 jgs 153 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 jgs 110
389 jgs 153 /* 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 jgs 150 }
439 jgs 153 }
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 jgs 150 }
508 jgs 153 }
509     fprintf(fileHandle_p, "</CellData>\n");
510 jgs 110 }
511 jgs 153 /* 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 jgs 113 /* finish off the piece */
642     fprintf(fileHandle_p, "</Piece>\n");
643 jgs 110
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