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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1028 - (hide annotations)
Wed Mar 14 00:15:24 2007 UTC (12 years, 7 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 jgs 150 /*
2 elspeth 616 ************************************************************
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 jgs 150 */
12 jgs 82
13 jgs 150
14 jgs 82 /**************************************************************/
15    
16     /* writes data and mesh in an opendx file */
17    
18     /**************************************************************/
19    
20     /* Author: Lutz Gross, gross@access.edu.au */
21 jgs 150 /* Version: $Id$ */
22 jgs 82
23     /**************************************************************/
24    
25     #include "Mesh.h"
26    
27 jgs 150 /**************************************************************/
28    
29 jgs 153 void Finley_Mesh_saveDX(const char * filename_p, Finley_Mesh *mesh_p, const dim_t num_data,char* *names_p,escriptDataC* *data_pp) {
30 gross 1028 char error_msg[LenErrorMsg_MAX], elemTypeStr[32];
31 jgs 82 /* 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 gross 1028 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 jgs 82 /* open the file and check handel */
47 jgs 153
48 gross 1028 /* 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 jgs 82 if (fileHandle_p==NULL) {
55 jgs 153 sprintf(error_msg,"File %s could not be opened for writing.",filename_p);
56 gross 1028 MEMFREE(isCellCentered);
57     fclose(fileHandle_p);
58 jgs 153 Finley_setError(IO_ERROR,error_msg);
59     return;
60 jgs 82 }
61 jgs 153 /* find the mesh type to be written */
62 gross 1028 nodetype=FINLEY_DEGREES_OF_FREEDOM;
63     elementtype=FINLEY_UNKNOWN;
64 jgs 153 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 gross 1028 MEMFREE(isCellCentered);
74     fclose(fileHandle_p);
75 jgs 153 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 gross 1028 MEMFREE(isCellCentered);
86     fclose(fileHandle_p);
87 jgs 153 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 gross 1028 MEMFREE(isCellCentered);
98     fclose(fileHandle_p);
99 jgs 153 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 gross 1028 MEMFREE(isCellCentered);
110     fclose(fileHandle_p);
111 jgs 153 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 gross 1028 MEMFREE(isCellCentered);
122     fclose(fileHandle_p);
123 jgs 153 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 gross 1028 MEMFREE(isCellCentered);
134     fclose(fileHandle_p);
135 jgs 153 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 gross 1028 MEMFREE(isCellCentered);
146     fclose(fileHandle_p);
147 jgs 153 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 gross 1028 MEMFREE(isCellCentered);
158     fclose(fileHandle_p);
159 jgs 153 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 gross 1028 MEMFREE(isCellCentered);
167     fclose(fileHandle_p);
168 jgs 153 return;
169     }
170     }
171     }
172     /* select number of points and the mesh component */
173 gross 1028 numPoints = mesh_p->Nodes->numNodes;
174     nDim = mesh_p->Nodes->numDim;
175 jgs 153 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 gross 1028 MEMFREE(isCellCentered);
198     fclose(fileHandle_p);
199 jgs 153 return;
200     }
201    
202     /* map finley element type to DX element type */
203 gross 1028 TypeId = elements->ReferenceElement->Type->TypeId;
204     numDXNodesPerElement=0;
205     numCells = elements->numElements;
206 jgs 153 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 gross 1028 MEMFREE(isCellCentered);
230     fclose(fileHandle_p);
231 jgs 153 return;
232     }
233    
234 jgs 82 /* positions */
235 jgs 153 fprintf(fileHandle_p, "object 1 class array type float rank 1 shape %d items %d data follows\n",nDim, mesh_p->Nodes->reducedNumNodes);
236 jgs 82 for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
237     if (mesh_p->Nodes->toReduced[i]>=0) {
238 jgs 153 for (j = 0; j < nDim; j++) fprintf(fileHandle_p, " %g",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);
239 jgs 82 fprintf(fileHandle_p, "\n");
240     }
241     }
242 jgs 153 /* connection table */
243 gross 1028 NN=elements->ReferenceElement->Type->numNodes;
244 jgs 153 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 jgs 82
252     /* data */
253 gross 1028 object_count=2;
254 jgs 153 for (i_data =0 ;i_data<num_data;++i_data) {
255     if (! isEmpty(data_pp[i_data])) {
256     object_count++;
257 gross 1028 rank=getDataPointRank(data_pp[i_data]);
258     nComp=getDataPointSize(data_pp[i_data]);
259 jgs 153 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 gross 903 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 jgs 153 }
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 jgs 82 }
301     }
302 jgs 153 fprintf(fileHandle_p, "attribute \"dep\" string \"positions\"\n");
303 jgs 82 }
304 jgs 153 }
305 jgs 82 }
306    
307     /* and finish it up */
308 jgs 153 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 jgs 82 fprintf(fileHandle_p, "end\n");
325     /* close the file */
326     fclose(fileHandle_p);
327 gross 1028 MEMFREE(isCellCentered);
328 jgs 82 return;
329     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26