/[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 471 - (hide annotations)
Fri Jan 27 01:33:02 2006 UTC (13 years, 9 months ago) by jgs
File MIME type: text/plain
File size: 13202 byte(s)
reorganise finley src tree to remove inc dir and src/finley directory

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26