/[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 110 - (hide annotations)
Mon Feb 14 04:14:42 2005 UTC (14 years, 6 months ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Mesh_saveVTK.c
File MIME type: text/plain
File size: 10028 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    
19     void Finley_Mesh_saveVTK(const char * filename_p, Finley_Mesh *mesh_p, escriptDataC* data_p) {
20     /* if there is no mesh we just return */
21     if (mesh_p==NULL) return;
22     /* some tables needed for reordering */
23     int resort[6][9]={
24     {0,1}, /* line */
25     {0,1,2}, /* triangle */
26     {0,1,2,3}, /* tetrahedron */
27     {0,3,1,2}, /* quadrilateral */
28     {3,0,7,4,2,1,6,5}, /* hexahedron */
29     };
30     Finley_ElementFile* elements=NULL;
31     char elemTypeStr[32];
32     int i,j,k,numVTKNodesPerElement,*resortIndex,isCellCentered,nodetype;
33     double* values,rtmp;
34     int nDim = mesh_p->Nodes->numDim;
35     /* open the file and check handle */
36     FILE * fileHandle_p = fopen(filename_p, "w");
37     if (fileHandle_p==NULL) {
38     Finley_ErrorCode=IO_ERROR;
39     sprintf(Finley_ErrorMsg,"File %s could not be opened for writing.",filename_p);
40     return;
41     }
42     /* xml header */
43     fprintf(fileHandle_p, "<?xml version=\"1.0\"?>\n");
44     fprintf(fileHandle_p,
45     "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\">\n");
46    
47     /* finley uses an unstructured mesh, so UnstructuredGrid *should* work */
48     fprintf(fileHandle_p, "<UnstructuredGrid>\n");
49    
50     /* is there only one "piece" to the data?? */
51     /* fprintf(fileHandle_p,
52     "<Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n",);
53     */
54    
55     /* now for the point data */
56     /* fprintf(fileHandle_p,
57     "<PointData Scalars=\"%s\" Vectors=\"%s\">\n",);
58     */
59    
60     fprintf(fileHandle_p, "</PointData>\n");
61    
62     /* now for the cell data */
63     /* fprintf(fileHandle_p,
64     "<CellData Scalars=\"%s\" Vectors=\"%s\">\n",);
65     */
66    
67     fprintf(fileHandle_p, "</CellData>\n");
68    
69     /* now for the points */
70     fprintf(fileHandle_p, "<Points>\n");
71     fprintf(fileHandle_p, "<DataArray NumberOfComponents=\"%d\">\n");
72     fprintf(fileHandle_p, "</DataArray>\n");
73     fprintf(fileHandle_p, "</Points>\n");
74    
75     /* now for the cells */
76     fprintf(fileHandle_p, "<Cells>\n");
77     /* fprintf(fileHandle_p,
78     "<DataArray type=\"%s\" Name=\"%s\" format=\"ascii\">\n",);
79     */
80     fprintf(fileHandle_p,
81     "</DataArray>\n");
82     fprintf(fileHandle_p, "</Cells>\n");
83    
84     /* finish off the piece */
85     fprintf(fileHandle_p, "</Piece>\n");
86    
87     /* positions */
88     fprintf(fileHandle_p, "object 1 class array type float rank 1 shape %d items %d data follows\n", nDim, mesh_p->Nodes->reducedNumNodes);
89     for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
90     if (mesh_p->Nodes->toReduced[i]>=0) {
91     fprintf(fileHandle_p, "%e", mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)]);
92     for (j = 1; j < nDim; j++) fprintf(fileHandle_p, " %f",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);
93     fprintf(fileHandle_p, "\n");
94     }
95     }
96     /* connections */
97     /* get a pointer to the relevant mesh component */
98     if (isEmpty(data_p)) {
99     elements=mesh_p->Elements;
100     } else {
101     switch(getFunctionSpaceType(data_p)) {
102     case(FINLEY_DEGREES_OF_FREEDOM):
103     nodetype=FINLEY_DEGREES_OF_FREEDOM;
104     case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):
105     nodetype=FINLEY_REDUCED_DEGREES_OF_FREEDOM;
106     case(FINLEY_NODES):
107     nodetype=FINLEY_NODES;
108     isCellCentered=FALSE;
109     elements=mesh_p->Elements;
110     break;
111     case(FINLEY_ELEMENTS):
112     isCellCentered=TRUE;
113     elements=mesh_p->Elements;
114     break;
115     case(FINLEY_FACE_ELEMENTS):
116     isCellCentered=TRUE;
117     elements=mesh_p->FaceElements;
118     break;
119     case(FINLEY_POINTS):
120     isCellCentered=TRUE;
121     elements=mesh_p->Points;
122     break;
123     case(FINLEY_CONTACT_ELEMENTS_1):
124     case(FINLEY_CONTACT_ELEMENTS_2):
125     isCellCentered=TRUE;
126     elements=mesh_p->ContactElements;
127     break;
128     default:
129     Finley_ErrorCode=TYPE_ERROR;
130     sprintf(Finley_ErrorMsg,"Finley does not know anything about function space type %d",getFunctionSpaceType(data_p));
131     break;
132     }
133     }
134     /* if no element table is present jump over the connection table */
135     if (elements!=NULL) {
136     ElementTypeId TypeId = elements->ReferenceElement->Type->TypeId;
137     if (TypeId==Line2 || TypeId==Line3 || TypeId==Line4 ) {
138     numVTKNodesPerElement=2;
139     resortIndex=resort[0];
140     strcpy(elemTypeStr, "lines");
141     } else if (TypeId==Tri3 || TypeId==Tri6 || TypeId==Tri9 || TypeId==Tri10 ) {
142     numVTKNodesPerElement = 3;
143     resortIndex=resort[1];
144     strcpy(elemTypeStr, "triangles");
145     } else if (TypeId==Rec4 || TypeId==Rec8 || TypeId==Rec9 || TypeId==Rec12 || TypeId==Rec16 ) {
146     numVTKNodesPerElement = 4;
147     resortIndex=resort[3];
148     strcpy(elemTypeStr, "quads");
149     } else if (TypeId==Tet4 || TypeId==Tet10 || TypeId==Tet16 ) {
150     numVTKNodesPerElement = 4;
151     resortIndex=resort[2];
152     strcpy(elemTypeStr, "tetrahedra");
153     } else if (TypeId==Hex8 || TypeId==Hex20 || TypeId==Hex32 ) {
154     numVTKNodesPerElement = 8;
155     resortIndex=resort[4];
156     strcpy(elemTypeStr, "cubes");
157     } else {
158     Finley_ErrorCode=VALUE_ERROR;
159     sprintf(Finley_ErrorMsg, "Element type %s is not supported by VTK",elements->ReferenceElement->Type->Name);
160     return;
161     }
162     int NN=elements->ReferenceElement->Type->numNodes;
163     fprintf(fileHandle_p, "object 2 class array type int rank 1 shape %d items %d data follows\n",numVTKNodesPerElement, elements->numElements);
164     for (i = 0; i < elements->numElements; i++) {
165     fprintf(fileHandle_p,"%d",mesh_p->Nodes->toReduced[mesh_p->Elements->Nodes[INDEX2(resortIndex[0], i, NN)]]);
166     for (j = 1; j < numVTKNodesPerElement; j++) {
167     fprintf(fileHandle_p," %d",mesh_p->Nodes->toReduced[mesh_p->Elements->Nodes[INDEX2(resortIndex[j], i, NN)]]);
168     }
169     fprintf(fileHandle_p, "\n");
170     }
171     fprintf(fileHandle_p, "attribute \"element type\" string \"%s\"\n",elemTypeStr);
172     fprintf(fileHandle_p, "attribute \"ref\" string \"positions\"\n");
173    
174     }
175     /* data */
176     if (!isEmpty(data_p)) {
177     int rank=getDataPointRank(data_p);
178     int* shape=getDataPointShape(data_p);
179     int nComp=getDataPointSize(data_p);
180     fprintf(fileHandle_p, "object 3 class array type float rank %d ", rank);
181     if (0 < rank) {
182     fprintf(fileHandle_p, "shape ");
183     for (i = 0; i < rank; i++) fprintf(fileHandle_p, "%d ", shape[i]);
184     }
185     if (isCellCentered) {
186     int numPointsPerSample=elements->ReferenceElement->numQuadNodes;
187     if (numPointsPerSample) {
188     fprintf(fileHandle_p, "items %d data follows\n", elements->numElements);
189     for (i=0;i<elements->numElements;i++) {
190     values=getSampleData(data_p,i);
191     for (k=0;k<nComp;k++) {
192     rtmp=0.;
193     for (j=0;j<numPointsPerSample;j++) rtmp+=values[INDEX2(k,j,nComp)];
194     fprintf(fileHandle_p, " %f", rtmp/numPointsPerSample);
195     }
196     fprintf(fileHandle_p, "\n");
197     }
198     fprintf(fileHandle_p, "attribute \"dep\" string \"connections\"\n");
199     }
200     } else {
201     fprintf(fileHandle_p, "items %d data follows\n", mesh_p->Nodes->reducedNumNodes);
202     for (i=0;i<mesh_p->Nodes->numNodes;i++) {
203     if (mesh_p->Nodes->toReduced[i]>=0) {
204     switch (nodetype) {
205     case(FINLEY_DEGREES_OF_FREEDOM):
206     values=getSampleData(data_p,mesh_p->Nodes->degreeOfFreedom[i]);
207     break;
208     case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):
209     values=getSampleData(data_p,mesh_p->Nodes->reducedDegreeOfFreedom[i]);
210     break;
211     case(FINLEY_NODES):
212     values=getSampleData(data_p,i);
213     break;
214     }
215     }
216     for (k=0;k<nComp;k++) fprintf(fileHandle_p, " %f", values[k]);
217     fprintf(fileHandle_p, "\n");
218     }
219     fprintf(fileHandle_p, "attribute \"dep\" string \"positions\"\n");
220     }
221     }
222    
223     /* and finish it up */
224     fprintf(fileHandle_p, "object 4 class field\n");
225     fprintf(fileHandle_p, "component \"positions\" value 1\n");
226     if (elements!=NULL) fprintf(fileHandle_p, "component \"connections\" value 2\n");
227     if (!isEmpty(data_p)) fprintf(fileHandle_p, "component \"data\" value 3\n");
228     fprintf(fileHandle_p, "end\n");
229    
230     fprintf(fileHandle_p, "</UnstructuredGrid>\n");
231     /* write the xml footer */
232     fprintf(fileHandle_p, "</VTKFile>\n");
233     /* close the file */
234     fclose(fileHandle_p);
235     return;
236     }
237    
238     /*
239     * $Log$
240     * Revision 1.2 2005/02/14 04:14:42 jgs
241     * *** empty log message ***
242     *
243     * Revision 1.1.2.2 2005/02/10 01:34:22 cochrane
244     * Quick fix to make sure that saveVTK compiles so that finley is still buildable. Apologies to those this has affected.
245     *
246     * Revision 1.1.2.1 2005/02/09 06:53:15 cochrane
247     * Initial import to repository. This is the file to implement saving finley/escript meshes out to vtk formatted files. It is basically just a hack of the opendx equivalent, with a lot of the opendx stuff still in the file, so it doesn't actually work just yet, but it probably needs to be added to the cvs.
248     *
249     * Revision 1.1.1.1 2004/10/26 06:53:57 jgs
250     * initial import of project esys2
251     *
252     * Revision 1.1 2004/07/27 08:27:11 gross
253     * Finley: saveDX added: now it is possible to write data on boundary and contact elements
254     *
255     */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26