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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 903 - (hide annotations)
Fri Nov 17 01:59:49 2006 UTC (13 years ago) by gross
Original Path: trunk/finley/src/Mesh_saveDX.c
File MIME type: text/plain
File size: 13081 byte(s)
bug with tagged data in vtk and dx writer fixed
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 jgs 150 char error_msg[LenErrorMsg_MAX];
31 jgs 82 /* if there is no mesh we just return */
32     if (mesh_p==NULL) return;
33 jgs 153
34 jgs 82 /* 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 jgs 153 int i,j,k,i_data;
43 jgs 82 /* open the file and check handel */
44 jgs 153
45 jgs 82 FILE * fileHandle_p = fopen(filename_p, "w");
46     if (fileHandle_p==NULL) {
47 jgs 153 sprintf(error_msg,"File %s could not be opened for writing.",filename_p);
48     Finley_setError(IO_ERROR,error_msg);
49     return;
50 jgs 82 }
51 jgs 153 /* 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 jgs 82 /* positions */
207 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);
208 jgs 82 for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
209     if (mesh_p->Nodes->toReduced[i]>=0) {
210 jgs 153 for (j = 0; j < nDim; j++) fprintf(fileHandle_p, " %g",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);
211 jgs 82 fprintf(fileHandle_p, "\n");
212     }
213     }
214 jgs 153 /* 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 jgs 82
224     /* data */
225 jgs 153 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 gross 903 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 jgs 153 }
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 jgs 82 }
274     }
275 jgs 153 fprintf(fileHandle_p, "attribute \"dep\" string \"positions\"\n");
276 jgs 82 }
277 jgs 153 }
278 jgs 82 }
279    
280     /* and finish it up */
281 jgs 153 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 jgs 82 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