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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 113 - (hide annotations)
Mon Feb 28 07:06:33 2005 UTC (14 years, 4 months ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Mesh_saveVTK.c
File MIME type: text/plain
File size: 17138 byte(s)
*** empty log message ***

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26