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

Contents of /trunk/finley/src/Mesh_saveDX.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 903 - (show annotations)
Fri Nov 17 01:59:49 2006 UTC (12 years, 11 months ago) by gross
File MIME type: text/plain
File size: 13081 byte(s)
bug with tagged data in vtk and dx writer fixed
1 /*
2 ************************************************************
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 */
12
13
14 /**************************************************************/
15
16 /* writes data and mesh in an opendx file */
17
18 /**************************************************************/
19
20 /* Author: Lutz Gross, gross@access.edu.au */
21 /* Version: $Id$ */
22
23 /**************************************************************/
24
25 #include "Mesh.h"
26
27 /**************************************************************/
28
29 void Finley_Mesh_saveDX(const char * filename_p, Finley_Mesh *mesh_p, const dim_t num_data,char* *names_p,escriptDataC* *data_pp) {
30 char error_msg[LenErrorMsg_MAX];
31 /* if there is no mesh we just return */
32 if (mesh_p==NULL) return;
33
34 /* some tables needed for reordering */
35 int resort[6][9]={
36 {0,1}, /* line */
37 {0,1,2}, /* triagle */
38 {0,1,2,3}, /* tetrahedron */
39 {0,3,1,2}, /* quadrilateral */
40 {3,0,7,4,2,1,6,5}, /* hexahedron */
41 };
42 int i,j,k,i_data;
43 /* open the file and check handel */
44
45 FILE * fileHandle_p = fopen(filename_p, "w");
46 if (fileHandle_p==NULL) {
47 sprintf(error_msg,"File %s could not be opened for writing.",filename_p);
48 Finley_setError(IO_ERROR,error_msg);
49 return;
50 }
51 /* 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];
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,"saveDX: cannot write given data in single file.");
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,"saveDX: cannot write given data in single file.");
74 return;
75 }
76 isCellCentered[i_data]=FALSE;
77 break;
78 case FINLEY_NODES:
79 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
80 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
81 elementtype=FINLEY_ELEMENTS;
82 } else {
83 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");
84 return;
85 }
86 isCellCentered[i_data]=FALSE;
87 break;
88 case FINLEY_ELEMENTS:
89 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
90 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
91 elementtype=FINLEY_ELEMENTS;
92 } else {
93 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");
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,"saveDX: cannot write given data in single file.");
104 return;
105 }
106 isCellCentered[i_data]=TRUE;
107 break;
108 case FINLEY_POINTS:
109 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
110 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_POINTS) {
111 elementtype=FINLEY_POINTS;
112 } else {
113 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");
114 return;
115 }
116 isCellCentered[i_data]=TRUE;
117 break;
118 case FINLEY_CONTACT_ELEMENTS_1:
119 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
120 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {
121 elementtype=FINLEY_CONTACT_ELEMENTS_1;
122 } else {
123 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");
124 return;
125 }
126 isCellCentered[i_data]=TRUE;
127 break;
128 case FINLEY_CONTACT_ELEMENTS_2:
129 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
130 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {
131 elementtype=FINLEY_CONTACT_ELEMENTS_1;
132 } else {
133 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");
134 return;
135 }
136 isCellCentered[i_data]=TRUE;
137 break;
138 default:
139 sprintf(error_msg,"saveDX: Finley does not know anything about function space type %d",getFunctionSpaceType(data_pp[i_data]));
140 Finley_setError(TYPE_ERROR,error_msg);
141 return;
142 }
143 }
144 }
145 /* select number of points and the mesh component */
146 int numPoints = mesh_p->Nodes->numNodes;
147 int nDim = mesh_p->Nodes->numDim;
148 if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
149 numPoints = mesh_p->Nodes->reducedNumNodes;
150 } else {
151 numPoints = mesh_p->Nodes->numNodes;
152 }
153 if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;
154 Finley_ElementFile* elements=NULL;
155 switch(elementtype) {
156 case FINLEY_ELEMENTS:
157 elements=mesh_p->Elements;
158 break;
159 case FINLEY_FACE_ELEMENTS:
160 elements=mesh_p->FaceElements;
161 break;
162 case FINLEY_POINTS:
163 elements=mesh_p->Points;
164 break;
165 case FINLEY_CONTACT_ELEMENTS_1:
166 elements=mesh_p->ContactElements;
167 break;
168 }
169 if (elements==NULL) {
170 Finley_setError(SYSTEM_ERROR,"saveDX: undefined element file");
171 return;
172 }
173
174 /* map finley element type to DX element type */
175 ElementTypeId TypeId = elements->ReferenceElement->Type->TypeId;
176 int *resortIndex=NULL;
177 int numDXNodesPerElement=0;
178 int numCells = elements->numElements;
179 char elemTypeStr[32];
180 if (TypeId==Line2 || TypeId==Line3 || TypeId==Line4 ) {
181 numDXNodesPerElement=2;
182 resortIndex=resort[0];
183 strcpy(elemTypeStr, "lines");
184 } else if (TypeId==Tri3 || TypeId==Tri6 || TypeId==Tri9 || TypeId==Tri10 ) {
185 numDXNodesPerElement = 3;
186 resortIndex=resort[1];
187 strcpy(elemTypeStr, "triangles");
188 } else if (TypeId==Rec4 || TypeId==Rec8 || TypeId==Rec9 || TypeId==Rec12 || TypeId==Rec16 ) {
189 numDXNodesPerElement = 4;
190 resortIndex=resort[3];
191 strcpy(elemTypeStr, "quads");
192 } else if (TypeId==Tet4 || TypeId==Tet10 || TypeId==Tet16 ) {
193 numDXNodesPerElement = 4;
194 resortIndex=resort[2];
195 strcpy(elemTypeStr, "tetrahedra");
196 } else if (TypeId==Hex8 || TypeId==Hex20 || TypeId==Hex32 ) {
197 numDXNodesPerElement = 8;
198 resortIndex=resort[4];
199 strcpy(elemTypeStr, "cubes");
200 } else {
201 sprintf(error_msg,"saveDX: Element type %s is not supported by DX",elements->ReferenceElement->Type->Name);
202 Finley_setError(VALUE_ERROR,error_msg);
203 return;
204 }
205
206 /* positions */
207 fprintf(fileHandle_p, "object 1 class array type float rank 1 shape %d items %d data follows\n",nDim, mesh_p->Nodes->reducedNumNodes);
208 for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
209 if (mesh_p->Nodes->toReduced[i]>=0) {
210 for (j = 0; j < nDim; j++) fprintf(fileHandle_p, " %g",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);
211 fprintf(fileHandle_p, "\n");
212 }
213 }
214 /* connection table */
215 int NN=elements->ReferenceElement->Type->numNodes;
216 fprintf(fileHandle_p, "object 2 class array type int rank 1 shape %d items %d data follows\n",numDXNodesPerElement, numCells);
217 for (i = 0; i < numCells; i++) {
218 for (j = 0; j < numDXNodesPerElement; j++) fprintf(fileHandle_p," %d",mesh_p->Nodes->toReduced[elements->Nodes[INDEX2(resortIndex[j], i, NN)]]);
219 fprintf(fileHandle_p, "\n");
220 }
221 fprintf(fileHandle_p, "attribute \"element type\" string \"%s\"\n",elemTypeStr);
222 fprintf(fileHandle_p, "attribute \"ref\" string \"positions\"\n");
223
224 /* data */
225 int object_count=2;
226 for (i_data =0 ;i_data<num_data;++i_data) {
227 if (! isEmpty(data_pp[i_data])) {
228 object_count++;
229 int rank=getDataPointRank(data_pp[i_data]);
230 int nComp=getDataPointSize(data_pp[i_data]);
231 double* values,rtmp;
232 fprintf(fileHandle_p, "object %d class array type float rank %d ",object_count,rank);
233 if (0 < rank) {
234 fprintf(fileHandle_p, "shape ");
235 for (i = 0; i < rank; i++) fprintf(fileHandle_p, "%d ", getDataPointShape(data_pp[i_data],i));
236 }
237 if (isCellCentered[i_data]) {
238 int numPointsPerSample=elements->ReferenceElement->numQuadNodes;
239 if (numPointsPerSample>0) {
240 fprintf(fileHandle_p, "items %d data follows\n", numCells);
241 for (i=0;i<elements->numElements;i++) {
242 values=getSampleData(data_pp[i_data],i);
243 for (k=0;k<nComp;k++) {
244 if ( isExpanded(data_pp[i_data]) ) {
245 rtmp=0.;
246 for (j=0;j<numPointsPerSample;j++) rtmp+=values[INDEX2(k,j,nComp)];
247 fprintf(fileHandle_p, " %g", rtmp/numPointsPerSample);
248 } else {
249 fprintf(fileHandle_p, " %g", values[k]);
250 }
251 }
252 fprintf(fileHandle_p, "\n");
253 }
254 fprintf(fileHandle_p, "attribute \"dep\" string \"connections\"\n");
255 }
256 } else {
257 fprintf(fileHandle_p, "items %d data follows\n", mesh_p->Nodes->reducedNumNodes);
258 for (i=0;i<mesh_p->Nodes->numNodes;i++) {
259 if (mesh_p->Nodes->toReduced[i]>=0) {
260 switch (getFunctionSpaceType(data_pp[i_data])) {
261 case FINLEY_DEGREES_OF_FREEDOM:
262 values=getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);
263 break;
264 case FINLEY_REDUCED_DEGREES_OF_FREEDOM:
265 values=getSampleData(data_pp[i_data],mesh_p->Nodes->reducedDegreeOfFreedom[i]);
266 break;
267 case FINLEY_NODES:
268 values=getSampleData(data_pp[i_data],i);
269 break;
270 }
271 for (k=0;k<nComp;k++) fprintf(fileHandle_p, " %g", values[k]);
272 fprintf(fileHandle_p, "\n");
273 }
274 }
275 fprintf(fileHandle_p, "attribute \"dep\" string \"positions\"\n");
276 }
277 }
278 }
279
280 /* and finish it up */
281 if (num_data==0) {
282 fprintf(fileHandle_p, "object %d class field\n",object_count+1);
283 fprintf(fileHandle_p, "component \"positions\" value 1\n");
284 fprintf(fileHandle_p, "component \"connections\" value 2\n");
285 } else {
286 object_count=2;
287 for (i_data=0; i_data<num_data;++i_data) {
288 if (! isEmpty(data_pp[i_data])) {
289 object_count++;
290 fprintf(fileHandle_p, "object \"%s\" class field\n",names_p[i_data]);
291 fprintf(fileHandle_p, "component \"positions\" value 1\n");
292 fprintf(fileHandle_p, "component \"connections\" value 2\n");
293 fprintf(fileHandle_p, "component \"data\" value %d\n",object_count);
294 }
295 }
296 }
297 fprintf(fileHandle_p, "end\n");
298 /* close the file */
299 fclose(fileHandle_p);
300 return;
301 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26