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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26