/[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 1028 - (show annotations)
Wed Mar 14 00:15:24 2007 UTC (12 years, 6 months ago) by gross
File MIME type: text/plain
File size: 14175 byte(s)
modifications to be compliant with _WIN32. The substitutes for asinh, acosh, atanh are still missing (erf will through an exception)
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], elemTypeStr[32];
31 /* some tables needed for reordering */
32 int resort[6][9]={
33 {0,1}, /* line */
34 {0,1,2}, /* triagle */
35 {0,1,2,3}, /* tetrahedron */
36 {0,3,1,2}, /* quadrilateral */
37 {3,0,7,4,2,1,6,5}, /* hexahedron */
38 };
39 FILE * fileHandle_p = NULL;
40 int i,j,k,i_data, nodetype, elementtype, numPoints = 0, nDim, *resortIndex=NULL,
41 numDXNodesPerElement=0, numCells, NN, object_count, rank, nComp;
42 double* values,rtmp;
43 bool_t *isCellCentered=NULL;
44 Finley_ElementFile* elements=NULL;
45 ElementTypeId TypeId;
46 /* open the file and check handel */
47
48 /* if there is no mesh we just return */
49 if (mesh_p==NULL) return;
50 isCellCentered=MEMALLOC(num_data, bool_t);
51 if (Finley_checkPtr(isCellCentered)) return;
52
53 fileHandle_p = fopen(filename_p, "w");
54 if (fileHandle_p==NULL) {
55 sprintf(error_msg,"File %s could not be opened for writing.",filename_p);
56 MEMFREE(isCellCentered);
57 fclose(fileHandle_p);
58 Finley_setError(IO_ERROR,error_msg);
59 return;
60 }
61 /* find the mesh type to be written */
62 nodetype=FINLEY_DEGREES_OF_FREEDOM;
63 elementtype=FINLEY_UNKNOWN;
64 for (i_data=0;i_data<num_data;++i_data) {
65 if (! isEmpty(data_pp[i_data])) {
66 switch(getFunctionSpaceType(data_pp[i_data])) {
67 case FINLEY_DEGREES_OF_FREEDOM:
68 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
69 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
70 elementtype=FINLEY_ELEMENTS;
71 } else {
72 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");
73 MEMFREE(isCellCentered);
74 fclose(fileHandle_p);
75 return;
76 }
77 isCellCentered[i_data]=FALSE;
78 break;
79 case FINLEY_REDUCED_DEGREES_OF_FREEDOM:
80 nodetype = FINLEY_REDUCED_DEGREES_OF_FREEDOM;
81 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
82 elementtype=FINLEY_ELEMENTS;
83 } else {
84 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");
85 MEMFREE(isCellCentered);
86 fclose(fileHandle_p);
87 return;
88 }
89 isCellCentered[i_data]=FALSE;
90 break;
91 case FINLEY_NODES:
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,"saveDX: cannot write given data in single file.");
97 MEMFREE(isCellCentered);
98 fclose(fileHandle_p);
99 return;
100 }
101 isCellCentered[i_data]=FALSE;
102 break;
103 case FINLEY_ELEMENTS:
104 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
105 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
106 elementtype=FINLEY_ELEMENTS;
107 } else {
108 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");
109 MEMFREE(isCellCentered);
110 fclose(fileHandle_p);
111 return;
112 }
113 isCellCentered[i_data]=TRUE;
114 break;
115 case FINLEY_FACE_ELEMENTS:
116 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
117 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_FACE_ELEMENTS) {
118 elementtype=FINLEY_FACE_ELEMENTS;
119 } else {
120 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");
121 MEMFREE(isCellCentered);
122 fclose(fileHandle_p);
123 return;
124 }
125 isCellCentered[i_data]=TRUE;
126 break;
127 case FINLEY_POINTS:
128 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
129 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_POINTS) {
130 elementtype=FINLEY_POINTS;
131 } else {
132 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");
133 MEMFREE(isCellCentered);
134 fclose(fileHandle_p);
135 return;
136 }
137 isCellCentered[i_data]=TRUE;
138 break;
139 case FINLEY_CONTACT_ELEMENTS_1:
140 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
141 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {
142 elementtype=FINLEY_CONTACT_ELEMENTS_1;
143 } else {
144 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");
145 MEMFREE(isCellCentered);
146 fclose(fileHandle_p);
147 return;
148 }
149 isCellCentered[i_data]=TRUE;
150 break;
151 case FINLEY_CONTACT_ELEMENTS_2:
152 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
153 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {
154 elementtype=FINLEY_CONTACT_ELEMENTS_1;
155 } else {
156 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");
157 MEMFREE(isCellCentered);
158 fclose(fileHandle_p);
159 return;
160 }
161 isCellCentered[i_data]=TRUE;
162 break;
163 default:
164 sprintf(error_msg,"saveDX: Finley does not know anything about function space type %d",getFunctionSpaceType(data_pp[i_data]));
165 Finley_setError(TYPE_ERROR,error_msg);
166 MEMFREE(isCellCentered);
167 fclose(fileHandle_p);
168 return;
169 }
170 }
171 }
172 /* select number of points and the mesh component */
173 numPoints = mesh_p->Nodes->numNodes;
174 nDim = mesh_p->Nodes->numDim;
175 if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
176 numPoints = mesh_p->Nodes->reducedNumNodes;
177 } else {
178 numPoints = mesh_p->Nodes->numNodes;
179 }
180 if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;
181 switch(elementtype) {
182 case FINLEY_ELEMENTS:
183 elements=mesh_p->Elements;
184 break;
185 case FINLEY_FACE_ELEMENTS:
186 elements=mesh_p->FaceElements;
187 break;
188 case FINLEY_POINTS:
189 elements=mesh_p->Points;
190 break;
191 case FINLEY_CONTACT_ELEMENTS_1:
192 elements=mesh_p->ContactElements;
193 break;
194 }
195 if (elements==NULL) {
196 Finley_setError(SYSTEM_ERROR,"saveDX: undefined element file");
197 MEMFREE(isCellCentered);
198 fclose(fileHandle_p);
199 return;
200 }
201
202 /* map finley element type to DX element type */
203 TypeId = elements->ReferenceElement->Type->TypeId;
204 numDXNodesPerElement=0;
205 numCells = elements->numElements;
206 if (TypeId==Line2 || TypeId==Line3 || TypeId==Line4 ) {
207 numDXNodesPerElement=2;
208 resortIndex=resort[0];
209 strcpy(elemTypeStr, "lines");
210 } else if (TypeId==Tri3 || TypeId==Tri6 || TypeId==Tri9 || TypeId==Tri10 ) {
211 numDXNodesPerElement = 3;
212 resortIndex=resort[1];
213 strcpy(elemTypeStr, "triangles");
214 } else if (TypeId==Rec4 || TypeId==Rec8 || TypeId==Rec9 || TypeId==Rec12 || TypeId==Rec16 ) {
215 numDXNodesPerElement = 4;
216 resortIndex=resort[3];
217 strcpy(elemTypeStr, "quads");
218 } else if (TypeId==Tet4 || TypeId==Tet10 || TypeId==Tet16 ) {
219 numDXNodesPerElement = 4;
220 resortIndex=resort[2];
221 strcpy(elemTypeStr, "tetrahedra");
222 } else if (TypeId==Hex8 || TypeId==Hex20 || TypeId==Hex32 ) {
223 numDXNodesPerElement = 8;
224 resortIndex=resort[4];
225 strcpy(elemTypeStr, "cubes");
226 } else {
227 sprintf(error_msg,"saveDX: Element type %s is not supported by DX",elements->ReferenceElement->Type->Name);
228 Finley_setError(VALUE_ERROR,error_msg);
229 MEMFREE(isCellCentered);
230 fclose(fileHandle_p);
231 return;
232 }
233
234 /* positions */
235 fprintf(fileHandle_p, "object 1 class array type float rank 1 shape %d items %d data follows\n",nDim, mesh_p->Nodes->reducedNumNodes);
236 for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
237 if (mesh_p->Nodes->toReduced[i]>=0) {
238 for (j = 0; j < nDim; j++) fprintf(fileHandle_p, " %g",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);
239 fprintf(fileHandle_p, "\n");
240 }
241 }
242 /* connection table */
243 NN=elements->ReferenceElement->Type->numNodes;
244 fprintf(fileHandle_p, "object 2 class array type int rank 1 shape %d items %d data follows\n",numDXNodesPerElement, numCells);
245 for (i = 0; i < numCells; i++) {
246 for (j = 0; j < numDXNodesPerElement; j++) fprintf(fileHandle_p," %d",mesh_p->Nodes->toReduced[elements->Nodes[INDEX2(resortIndex[j], i, NN)]]);
247 fprintf(fileHandle_p, "\n");
248 }
249 fprintf(fileHandle_p, "attribute \"element type\" string \"%s\"\n",elemTypeStr);
250 fprintf(fileHandle_p, "attribute \"ref\" string \"positions\"\n");
251
252 /* data */
253 object_count=2;
254 for (i_data =0 ;i_data<num_data;++i_data) {
255 if (! isEmpty(data_pp[i_data])) {
256 object_count++;
257 rank=getDataPointRank(data_pp[i_data]);
258 nComp=getDataPointSize(data_pp[i_data]);
259 fprintf(fileHandle_p, "object %d class array type float rank %d ",object_count,rank);
260 if (0 < rank) {
261 fprintf(fileHandle_p, "shape ");
262 for (i = 0; i < rank; i++) fprintf(fileHandle_p, "%d ", getDataPointShape(data_pp[i_data],i));
263 }
264 if (isCellCentered[i_data]) {
265 int numPointsPerSample=elements->ReferenceElement->numQuadNodes;
266 if (numPointsPerSample>0) {
267 fprintf(fileHandle_p, "items %d data follows\n", numCells);
268 for (i=0;i<elements->numElements;i++) {
269 values=getSampleData(data_pp[i_data],i);
270 for (k=0;k<nComp;k++) {
271 if ( isExpanded(data_pp[i_data]) ) {
272 rtmp=0.;
273 for (j=0;j<numPointsPerSample;j++) rtmp+=values[INDEX2(k,j,nComp)];
274 fprintf(fileHandle_p, " %g", rtmp/numPointsPerSample);
275 } else {
276 fprintf(fileHandle_p, " %g", values[k]);
277 }
278 }
279 fprintf(fileHandle_p, "\n");
280 }
281 fprintf(fileHandle_p, "attribute \"dep\" string \"connections\"\n");
282 }
283 } else {
284 fprintf(fileHandle_p, "items %d data follows\n", mesh_p->Nodes->reducedNumNodes);
285 for (i=0;i<mesh_p->Nodes->numNodes;i++) {
286 if (mesh_p->Nodes->toReduced[i]>=0) {
287 switch (getFunctionSpaceType(data_pp[i_data])) {
288 case FINLEY_DEGREES_OF_FREEDOM:
289 values=getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);
290 break;
291 case FINLEY_REDUCED_DEGREES_OF_FREEDOM:
292 values=getSampleData(data_pp[i_data],mesh_p->Nodes->reducedDegreeOfFreedom[i]);
293 break;
294 case FINLEY_NODES:
295 values=getSampleData(data_pp[i_data],i);
296 break;
297 }
298 for (k=0;k<nComp;k++) fprintf(fileHandle_p, " %g", values[k]);
299 fprintf(fileHandle_p, "\n");
300 }
301 }
302 fprintf(fileHandle_p, "attribute \"dep\" string \"positions\"\n");
303 }
304 }
305 }
306
307 /* and finish it up */
308 if (num_data==0) {
309 fprintf(fileHandle_p, "object %d class field\n",object_count+1);
310 fprintf(fileHandle_p, "component \"positions\" value 1\n");
311 fprintf(fileHandle_p, "component \"connections\" value 2\n");
312 } else {
313 object_count=2;
314 for (i_data=0; i_data<num_data;++i_data) {
315 if (! isEmpty(data_pp[i_data])) {
316 object_count++;
317 fprintf(fileHandle_p, "object \"%s\" class field\n",names_p[i_data]);
318 fprintf(fileHandle_p, "component \"positions\" value 1\n");
319 fprintf(fileHandle_p, "component \"connections\" value 2\n");
320 fprintf(fileHandle_p, "component \"data\" value %d\n",object_count);
321 }
322 }
323 }
324 fprintf(fileHandle_p, "end\n");
325 /* close the file */
326 fclose(fileHandle_p);
327 MEMFREE(isCellCentered);
328 return;
329 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26