/[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 1811 - (hide annotations)
Thu Sep 25 23:11:13 2008 UTC (10 years, 10 months ago) by ksteube
File MIME type: text/plain
File size: 11649 byte(s)
Copyright updated in all files

1 jgs 82
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4     * Copyright (c) 2003-2008 by University of Queensland
5     * Earth Systems Science Computational Center (ESSCC)
6     * http://www.uq.edu.au/esscc
7     *
8     * Primary Business: Queensland, Australia
9     * Licensed under the Open Software License version 3.0
10     * http://www.opensource.org/licenses/osl-3.0.php
11     *
12     *******************************************************/
13 jgs 82
14 ksteube 1811
15 jgs 82 /**************************************************************/
16    
17 ksteube 1312 /* writes data and mesh in an opendx file */
18     /* the input data needs to be cell centered or on reducedNodes */
19 jgs 82
20     /**************************************************************/
21    
22     #include "Mesh.h"
23 gross 1062 #include "Assemble.h"
24 jgs 82
25 jgs 150 /**************************************************************/
26    
27 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) {
28 gross 1028 char error_msg[LenErrorMsg_MAX], elemTypeStr[32];
29 jgs 82 /* some tables needed for reordering */
30     int resort[6][9]={
31     {0,1}, /* line */
32     {0,1,2}, /* triagle */
33     {0,1,2,3}, /* tetrahedron */
34     {0,3,1,2}, /* quadrilateral */
35     {3,0,7,4,2,1,6,5}, /* hexahedron */
36     };
37 gross 1028 FILE * fileHandle_p = NULL;
38 ksteube 1312 int i,j,k,i_data, elementtype, numPoints = 0, nDim, *resortIndex=NULL, p,
39 gross 1062 numDXNodesPerElement=0, numCells, NN, object_count, rank, nComp, numPointsPerSample;
40 gross 1028 double* values,rtmp;
41     bool_t *isCellCentered=NULL;
42     Finley_ElementFile* elements=NULL;
43     ElementTypeId TypeId;
44 jgs 82 /* open the file and check handel */
45 jgs 153
46 gross 1028 /* if there is no mesh we just return */
47     if (mesh_p==NULL) return;
48     isCellCentered=MEMALLOC(num_data, bool_t);
49     if (Finley_checkPtr(isCellCentered)) return;
50    
51     fileHandle_p = fopen(filename_p, "w");
52 jgs 82 if (fileHandle_p==NULL) {
53 jgs 153 sprintf(error_msg,"File %s could not be opened for writing.",filename_p);
54 gross 1028 MEMFREE(isCellCentered);
55     fclose(fileHandle_p);
56 jgs 153 Finley_setError(IO_ERROR,error_msg);
57     return;
58 jgs 82 }
59 jgs 153 /* find the mesh type to be written */
60 gross 1338 elementtype=FINLEY_UNKNOWN;
61 jgs 153 for (i_data=0;i_data<num_data;++i_data) {
62     if (! isEmpty(data_pp[i_data])) {
63     switch(getFunctionSpaceType(data_pp[i_data])) {
64 ksteube 1312 case FINLEY_REDUCED_NODES:
65 jgs 153 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
66     elementtype=FINLEY_ELEMENTS;
67     } else {
68 gross 1361 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data into single file.");
69 gross 1028 MEMFREE(isCellCentered);
70     fclose(fileHandle_p);
71 jgs 153 return;
72     }
73     isCellCentered[i_data]=FALSE;
74     break;
75     case FINLEY_ELEMENTS:
76 gross 1062 case FINLEY_REDUCED_ELEMENTS:
77 jgs 153 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
78     elementtype=FINLEY_ELEMENTS;
79     } else {
80 gross 1361 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data into single file.");
81 gross 1028 MEMFREE(isCellCentered);
82     fclose(fileHandle_p);
83 jgs 153 return;
84     }
85     isCellCentered[i_data]=TRUE;
86     break;
87     case FINLEY_FACE_ELEMENTS:
88 gross 1062 case FINLEY_REDUCED_FACE_ELEMENTS:
89 jgs 153 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_FACE_ELEMENTS) {
90     elementtype=FINLEY_FACE_ELEMENTS;
91     } else {
92 gross 1361 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data into single file.");
93 gross 1028 MEMFREE(isCellCentered);
94     fclose(fileHandle_p);
95 jgs 153 return;
96     }
97     isCellCentered[i_data]=TRUE;
98     break;
99     case FINLEY_POINTS:
100     if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_POINTS) {
101     elementtype=FINLEY_POINTS;
102     } else {
103 gross 1361 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data into single file.");
104 gross 1028 MEMFREE(isCellCentered);
105     fclose(fileHandle_p);
106 jgs 153 return;
107     }
108     isCellCentered[i_data]=TRUE;
109     break;
110     case FINLEY_CONTACT_ELEMENTS_1:
111 gross 1062 case FINLEY_REDUCED_CONTACT_ELEMENTS_1:
112 jgs 153 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {
113     elementtype=FINLEY_CONTACT_ELEMENTS_1;
114     } else {
115 gross 1361 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data into single file.");
116 gross 1028 MEMFREE(isCellCentered);
117     fclose(fileHandle_p);
118 jgs 153 return;
119     }
120     isCellCentered[i_data]=TRUE;
121     break;
122     case FINLEY_CONTACT_ELEMENTS_2:
123 gross 1062 case FINLEY_REDUCED_CONTACT_ELEMENTS_2:
124 jgs 153 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {
125     elementtype=FINLEY_CONTACT_ELEMENTS_1;
126     } else {
127 gross 1361 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data into single file.");
128 gross 1028 MEMFREE(isCellCentered);
129     fclose(fileHandle_p);
130 jgs 153 return;
131     }
132     isCellCentered[i_data]=TRUE;
133     break;
134     default:
135 gross 1361 sprintf(error_msg,"saveDX: I don't know how to handle function space type %d",getFunctionSpaceType(data_pp[i_data]));
136 jgs 153 Finley_setError(TYPE_ERROR,error_msg);
137 gross 1028 MEMFREE(isCellCentered);
138     fclose(fileHandle_p);
139 jgs 153 return;
140     }
141     }
142     }
143     /* select number of points and the mesh component */
144 ksteube 1312 numPoints = mesh_p->Nodes->reducedNodesMapping->numTargets;
145 gross 1028 nDim = mesh_p->Nodes->numDim;
146 jgs 153 if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;
147     switch(elementtype) {
148     case FINLEY_ELEMENTS:
149     elements=mesh_p->Elements;
150     break;
151     case FINLEY_FACE_ELEMENTS:
152     elements=mesh_p->FaceElements;
153     break;
154     case FINLEY_POINTS:
155     elements=mesh_p->Points;
156     break;
157 gross 1062 case FINLEY_CONTACT_ELEMENTS_2:
158     elements=mesh_p->ContactElements;
159     break;
160 jgs 153 case FINLEY_CONTACT_ELEMENTS_1:
161     elements=mesh_p->ContactElements;
162     break;
163     }
164     if (elements==NULL) {
165     Finley_setError(SYSTEM_ERROR,"saveDX: undefined element file");
166 gross 1028 MEMFREE(isCellCentered);
167     fclose(fileHandle_p);
168 jgs 153 return;
169     }
170    
171     /* map finley element type to DX element type */
172 gross 1028 TypeId = elements->ReferenceElement->Type->TypeId;
173     numDXNodesPerElement=0;
174     numCells = elements->numElements;
175 jgs 153 if (TypeId==Line2 || TypeId==Line3 || TypeId==Line4 ) {
176     numDXNodesPerElement=2;
177     resortIndex=resort[0];
178     strcpy(elemTypeStr, "lines");
179     } else if (TypeId==Tri3 || TypeId==Tri6 || TypeId==Tri9 || TypeId==Tri10 ) {
180     numDXNodesPerElement = 3;
181     resortIndex=resort[1];
182     strcpy(elemTypeStr, "triangles");
183     } else if (TypeId==Rec4 || TypeId==Rec8 || TypeId==Rec9 || TypeId==Rec12 || TypeId==Rec16 ) {
184     numDXNodesPerElement = 4;
185     resortIndex=resort[3];
186     strcpy(elemTypeStr, "quads");
187     } else if (TypeId==Tet4 || TypeId==Tet10 || TypeId==Tet16 ) {
188     numDXNodesPerElement = 4;
189     resortIndex=resort[2];
190     strcpy(elemTypeStr, "tetrahedra");
191     } else if (TypeId==Hex8 || TypeId==Hex20 || TypeId==Hex32 ) {
192     numDXNodesPerElement = 8;
193     resortIndex=resort[4];
194     strcpy(elemTypeStr, "cubes");
195     } else {
196     sprintf(error_msg,"saveDX: Element type %s is not supported by DX",elements->ReferenceElement->Type->Name);
197     Finley_setError(VALUE_ERROR,error_msg);
198 gross 1028 MEMFREE(isCellCentered);
199     fclose(fileHandle_p);
200 jgs 153 return;
201     }
202    
203 jgs 82 /* positions */
204 ksteube 1312 fprintf(fileHandle_p, "object 1 class array type float rank 1 shape %d items %d data follows\n",nDim, numPoints);
205     for (i = 0; i < numPoints; i++) {
206     p=mesh_p->Nodes->reducedNodesMapping->map[i];
207     for (j = 0; j < nDim; j++) fprintf(fileHandle_p, " %g",mesh_p->Nodes->Coordinates[INDEX2(j, p, nDim)]);
208 jgs 82 fprintf(fileHandle_p, "\n");
209     }
210 jgs 153 /* connection table */
211 ksteube 1312 NN=elements->numNodes;
212 jgs 153 fprintf(fileHandle_p, "object 2 class array type int rank 1 shape %d items %d data follows\n",numDXNodesPerElement, numCells);
213     for (i = 0; i < numCells; i++) {
214 ksteube 1312 for (j = 0; j < numDXNodesPerElement; j++)
215     fprintf(fileHandle_p," %d",mesh_p->Nodes->reducedNodesMapping->target[elements->Nodes[INDEX2(resortIndex[j], i, NN)]]);
216 jgs 153 fprintf(fileHandle_p, "\n");
217     }
218     fprintf(fileHandle_p, "attribute \"element type\" string \"%s\"\n",elemTypeStr);
219     fprintf(fileHandle_p, "attribute \"ref\" string \"positions\"\n");
220 jgs 82
221     /* data */
222 gross 1028 object_count=2;
223 jgs 153 for (i_data =0 ;i_data<num_data;++i_data) {
224     if (! isEmpty(data_pp[i_data])) {
225     object_count++;
226 gross 1028 rank=getDataPointRank(data_pp[i_data]);
227     nComp=getDataPointSize(data_pp[i_data]);
228 jgs 153 fprintf(fileHandle_p, "object %d class array type float rank %d ",object_count,rank);
229     if (0 < rank) {
230     fprintf(fileHandle_p, "shape ");
231     for (i = 0; i < rank; i++) fprintf(fileHandle_p, "%d ", getDataPointShape(data_pp[i_data],i));
232     }
233     if (isCellCentered[i_data]) {
234 gross 1062 if (Finley_Assemble_reducedIntegrationOrder(data_pp[i_data])) {
235     numPointsPerSample=elements->ReferenceElementReducedOrder->numQuadNodes;
236     } else {
237     numPointsPerSample=elements->ReferenceElement->numQuadNodes;
238     }
239 jgs 153 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 ksteube 1312 fprintf(fileHandle_p, "items %d data follows\n", numPoints);
258     for (i=0;i<numPoints;i++) {
259     values=getSampleData(data_pp[i_data],i);
260     for (k=0;k<nComp;k++) fprintf(fileHandle_p, " %g", values[k]);
261     fprintf(fileHandle_p, "\n");
262 jgs 82 }
263 jgs 153 fprintf(fileHandle_p, "attribute \"dep\" string \"positions\"\n");
264 jgs 82 }
265 jgs 153 }
266 jgs 82 }
267    
268     /* and finish it up */
269 jgs 153 if (num_data==0) {
270     fprintf(fileHandle_p, "object %d class field\n",object_count+1);
271     fprintf(fileHandle_p, "component \"positions\" value 1\n");
272     fprintf(fileHandle_p, "component \"connections\" value 2\n");
273     } else {
274     object_count=2;
275     for (i_data=0; i_data<num_data;++i_data) {
276     if (! isEmpty(data_pp[i_data])) {
277     object_count++;
278     fprintf(fileHandle_p, "object \"%s\" class field\n",names_p[i_data]);
279     fprintf(fileHandle_p, "component \"positions\" value 1\n");
280     fprintf(fileHandle_p, "component \"connections\" value 2\n");
281     fprintf(fileHandle_p, "component \"data\" value %d\n",object_count);
282     }
283     }
284     }
285 jgs 82 fprintf(fileHandle_p, "end\n");
286     /* close the file */
287     fclose(fileHandle_p);
288 gross 1028 MEMFREE(isCellCentered);
289 jgs 82 return;
290     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26