/[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 1384 - (hide annotations)
Fri Jan 11 02:29:38 2008 UTC (11 years, 11 months ago) by phornby
Original Path: temp_trunk_copy/finley/src/Mesh_saveDX.c
File MIME type: text/plain
File size: 11684 byte(s)
Make a temp copy of the trunk before checking in the windows changes


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26