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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
  ViewVC Help
Powered by ViewVC 1.1.26