/[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 1062 - (hide annotations)
Mon Mar 26 06:17:53 2007 UTC (12 years, 6 months ago) by gross
File MIME type: text/plain
File size: 14697 byte(s)
reduced integration schemes are implemented now for grad, integrate, etc. Tests still to be added.
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 gross 1062 #include "Assemble.h"
27 jgs 82
28 jgs 150 /**************************************************************/
29    
30 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) {
31 gross 1028 char error_msg[LenErrorMsg_MAX], elemTypeStr[32];
32 jgs 82 /* some tables needed for reordering */
33     int resort[6][9]={
34     {0,1}, /* line */
35     {0,1,2}, /* triagle */
36     {0,1,2,3}, /* tetrahedron */
37     {0,3,1,2}, /* quadrilateral */
38     {3,0,7,4,2,1,6,5}, /* hexahedron */
39     };
40 gross 1028 FILE * fileHandle_p = NULL;
41     int i,j,k,i_data, nodetype, elementtype, numPoints = 0, nDim, *resortIndex=NULL,
42 gross 1062 numDXNodesPerElement=0, numCells, NN, object_count, rank, nComp, numPointsPerSample;
43 gross 1028 double* values,rtmp;
44     bool_t *isCellCentered=NULL;
45     Finley_ElementFile* elements=NULL;
46     ElementTypeId TypeId;
47 jgs 82 /* open the file and check handel */
48 jgs 153
49 gross 1028 /* if there is no mesh we just return */
50     if (mesh_p==NULL) return;
51     isCellCentered=MEMALLOC(num_data, bool_t);
52     if (Finley_checkPtr(isCellCentered)) return;
53    
54     fileHandle_p = fopen(filename_p, "w");
55 jgs 82 if (fileHandle_p==NULL) {
56 jgs 153 sprintf(error_msg,"File %s could not be opened for writing.",filename_p);
57 gross 1028 MEMFREE(isCellCentered);
58     fclose(fileHandle_p);
59 jgs 153 Finley_setError(IO_ERROR,error_msg);
60     return;
61 jgs 82 }
62 jgs 153 /* find the mesh type to be written */
63 gross 1028 nodetype=FINLEY_DEGREES_OF_FREEDOM;
64     elementtype=FINLEY_UNKNOWN;
65 jgs 153 for (i_data=0;i_data<num_data;++i_data) {
66     if (! isEmpty(data_pp[i_data])) {
67     switch(getFunctionSpaceType(data_pp[i_data])) {
68     case FINLEY_DEGREES_OF_FREEDOM:
69     nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_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 gross 1028 MEMFREE(isCellCentered);
75     fclose(fileHandle_p);
76 jgs 153 return;
77     }
78     isCellCentered[i_data]=FALSE;
79     break;
80     case FINLEY_REDUCED_DEGREES_OF_FREEDOM:
81     nodetype = FINLEY_REDUCED_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 gross 1028 MEMFREE(isCellCentered);
87     fclose(fileHandle_p);
88 jgs 153 return;
89     }
90     isCellCentered[i_data]=FALSE;
91     break;
92     case FINLEY_NODES:
93     nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
94     if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
95     elementtype=FINLEY_ELEMENTS;
96     } else {
97     Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");
98 gross 1028 MEMFREE(isCellCentered);
99     fclose(fileHandle_p);
100 jgs 153 return;
101     }
102     isCellCentered[i_data]=FALSE;
103     break;
104     case FINLEY_ELEMENTS:
105 gross 1062 case FINLEY_REDUCED_ELEMENTS:
106 jgs 153 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
107     if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
108     elementtype=FINLEY_ELEMENTS;
109     } else {
110     Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");
111 gross 1028 MEMFREE(isCellCentered);
112     fclose(fileHandle_p);
113 jgs 153 return;
114     }
115     isCellCentered[i_data]=TRUE;
116     break;
117     case FINLEY_FACE_ELEMENTS:
118 gross 1062 case FINLEY_REDUCED_FACE_ELEMENTS:
119 jgs 153 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
120     if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_FACE_ELEMENTS) {
121     elementtype=FINLEY_FACE_ELEMENTS;
122     } else {
123     Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");
124 gross 1028 MEMFREE(isCellCentered);
125     fclose(fileHandle_p);
126 jgs 153 return;
127     }
128     isCellCentered[i_data]=TRUE;
129     break;
130     case FINLEY_POINTS:
131     nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
132     if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_POINTS) {
133     elementtype=FINLEY_POINTS;
134     } else {
135     Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");
136 gross 1028 MEMFREE(isCellCentered);
137     fclose(fileHandle_p);
138 jgs 153 return;
139     }
140     isCellCentered[i_data]=TRUE;
141     break;
142     case FINLEY_CONTACT_ELEMENTS_1:
143 gross 1062 case FINLEY_REDUCED_CONTACT_ELEMENTS_1:
144 jgs 153 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
145     if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {
146     elementtype=FINLEY_CONTACT_ELEMENTS_1;
147     } else {
148     Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");
149 gross 1028 MEMFREE(isCellCentered);
150     fclose(fileHandle_p);
151 jgs 153 return;
152     }
153     isCellCentered[i_data]=TRUE;
154     break;
155     case FINLEY_CONTACT_ELEMENTS_2:
156 gross 1062 case FINLEY_REDUCED_CONTACT_ELEMENTS_2:
157 jgs 153 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
158     if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {
159     elementtype=FINLEY_CONTACT_ELEMENTS_1;
160     } else {
161     Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");
162 gross 1028 MEMFREE(isCellCentered);
163     fclose(fileHandle_p);
164 jgs 153 return;
165     }
166     isCellCentered[i_data]=TRUE;
167     break;
168     default:
169     sprintf(error_msg,"saveDX: Finley does not know anything about function space type %d",getFunctionSpaceType(data_pp[i_data]));
170     Finley_setError(TYPE_ERROR,error_msg);
171 gross 1028 MEMFREE(isCellCentered);
172     fclose(fileHandle_p);
173 jgs 153 return;
174     }
175     }
176     }
177     /* select number of points and the mesh component */
178 gross 1028 numPoints = mesh_p->Nodes->numNodes;
179     nDim = mesh_p->Nodes->numDim;
180 jgs 153 if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
181     numPoints = mesh_p->Nodes->reducedNumNodes;
182     } else {
183     numPoints = mesh_p->Nodes->numNodes;
184     }
185     if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;
186     switch(elementtype) {
187     case FINLEY_ELEMENTS:
188     elements=mesh_p->Elements;
189     break;
190     case FINLEY_FACE_ELEMENTS:
191     elements=mesh_p->FaceElements;
192     break;
193     case FINLEY_POINTS:
194     elements=mesh_p->Points;
195     break;
196 gross 1062 case FINLEY_CONTACT_ELEMENTS_2:
197     elements=mesh_p->ContactElements;
198     break;
199 jgs 153 case FINLEY_CONTACT_ELEMENTS_1:
200     elements=mesh_p->ContactElements;
201     break;
202     }
203     if (elements==NULL) {
204     Finley_setError(SYSTEM_ERROR,"saveDX: undefined element file");
205 gross 1028 MEMFREE(isCellCentered);
206     fclose(fileHandle_p);
207 jgs 153 return;
208     }
209    
210     /* map finley element type to DX element type */
211 gross 1028 TypeId = elements->ReferenceElement->Type->TypeId;
212     numDXNodesPerElement=0;
213     numCells = elements->numElements;
214 jgs 153 if (TypeId==Line2 || TypeId==Line3 || TypeId==Line4 ) {
215     numDXNodesPerElement=2;
216     resortIndex=resort[0];
217     strcpy(elemTypeStr, "lines");
218     } else if (TypeId==Tri3 || TypeId==Tri6 || TypeId==Tri9 || TypeId==Tri10 ) {
219     numDXNodesPerElement = 3;
220     resortIndex=resort[1];
221     strcpy(elemTypeStr, "triangles");
222     } else if (TypeId==Rec4 || TypeId==Rec8 || TypeId==Rec9 || TypeId==Rec12 || TypeId==Rec16 ) {
223     numDXNodesPerElement = 4;
224     resortIndex=resort[3];
225     strcpy(elemTypeStr, "quads");
226     } else if (TypeId==Tet4 || TypeId==Tet10 || TypeId==Tet16 ) {
227     numDXNodesPerElement = 4;
228     resortIndex=resort[2];
229     strcpy(elemTypeStr, "tetrahedra");
230     } else if (TypeId==Hex8 || TypeId==Hex20 || TypeId==Hex32 ) {
231     numDXNodesPerElement = 8;
232     resortIndex=resort[4];
233     strcpy(elemTypeStr, "cubes");
234     } else {
235     sprintf(error_msg,"saveDX: Element type %s is not supported by DX",elements->ReferenceElement->Type->Name);
236     Finley_setError(VALUE_ERROR,error_msg);
237 gross 1028 MEMFREE(isCellCentered);
238     fclose(fileHandle_p);
239 jgs 153 return;
240     }
241    
242 jgs 82 /* positions */
243 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);
244 jgs 82 for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
245     if (mesh_p->Nodes->toReduced[i]>=0) {
246 jgs 153 for (j = 0; j < nDim; j++) fprintf(fileHandle_p, " %g",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);
247 jgs 82 fprintf(fileHandle_p, "\n");
248     }
249     }
250 jgs 153 /* connection table */
251 gross 1028 NN=elements->ReferenceElement->Type->numNodes;
252 jgs 153 fprintf(fileHandle_p, "object 2 class array type int rank 1 shape %d items %d data follows\n",numDXNodesPerElement, numCells);
253     for (i = 0; i < numCells; i++) {
254     for (j = 0; j < numDXNodesPerElement; j++) fprintf(fileHandle_p," %d",mesh_p->Nodes->toReduced[elements->Nodes[INDEX2(resortIndex[j], i, NN)]]);
255     fprintf(fileHandle_p, "\n");
256     }
257     fprintf(fileHandle_p, "attribute \"element type\" string \"%s\"\n",elemTypeStr);
258     fprintf(fileHandle_p, "attribute \"ref\" string \"positions\"\n");
259 jgs 82
260     /* data */
261 gross 1028 object_count=2;
262 jgs 153 for (i_data =0 ;i_data<num_data;++i_data) {
263     if (! isEmpty(data_pp[i_data])) {
264     object_count++;
265 gross 1028 rank=getDataPointRank(data_pp[i_data]);
266     nComp=getDataPointSize(data_pp[i_data]);
267 jgs 153 fprintf(fileHandle_p, "object %d class array type float rank %d ",object_count,rank);
268     if (0 < rank) {
269     fprintf(fileHandle_p, "shape ");
270     for (i = 0; i < rank; i++) fprintf(fileHandle_p, "%d ", getDataPointShape(data_pp[i_data],i));
271     }
272     if (isCellCentered[i_data]) {
273 gross 1062 if (Finley_Assemble_reducedIntegrationOrder(data_pp[i_data])) {
274     numPointsPerSample=elements->ReferenceElementReducedOrder->numQuadNodes;
275     } else {
276     numPointsPerSample=elements->ReferenceElement->numQuadNodes;
277     }
278 jgs 153 if (numPointsPerSample>0) {
279     fprintf(fileHandle_p, "items %d data follows\n", numCells);
280     for (i=0;i<elements->numElements;i++) {
281     values=getSampleData(data_pp[i_data],i);
282     for (k=0;k<nComp;k++) {
283 gross 903 if ( isExpanded(data_pp[i_data]) ) {
284     rtmp=0.;
285     for (j=0;j<numPointsPerSample;j++) rtmp+=values[INDEX2(k,j,nComp)];
286     fprintf(fileHandle_p, " %g", rtmp/numPointsPerSample);
287     } else {
288     fprintf(fileHandle_p, " %g", values[k]);
289     }
290 jgs 153 }
291     fprintf(fileHandle_p, "\n");
292     }
293     fprintf(fileHandle_p, "attribute \"dep\" string \"connections\"\n");
294     }
295     } else {
296     fprintf(fileHandle_p, "items %d data follows\n", mesh_p->Nodes->reducedNumNodes);
297     for (i=0;i<mesh_p->Nodes->numNodes;i++) {
298     if (mesh_p->Nodes->toReduced[i]>=0) {
299     switch (getFunctionSpaceType(data_pp[i_data])) {
300     case FINLEY_DEGREES_OF_FREEDOM:
301     values=getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);
302     break;
303     case FINLEY_REDUCED_DEGREES_OF_FREEDOM:
304     values=getSampleData(data_pp[i_data],mesh_p->Nodes->reducedDegreeOfFreedom[i]);
305     break;
306     case FINLEY_NODES:
307     values=getSampleData(data_pp[i_data],i);
308     break;
309     }
310     for (k=0;k<nComp;k++) fprintf(fileHandle_p, " %g", values[k]);
311     fprintf(fileHandle_p, "\n");
312 jgs 82 }
313     }
314 jgs 153 fprintf(fileHandle_p, "attribute \"dep\" string \"positions\"\n");
315 jgs 82 }
316 jgs 153 }
317 jgs 82 }
318    
319     /* and finish it up */
320 jgs 153 if (num_data==0) {
321     fprintf(fileHandle_p, "object %d class field\n",object_count+1);
322     fprintf(fileHandle_p, "component \"positions\" value 1\n");
323     fprintf(fileHandle_p, "component \"connections\" value 2\n");
324     } else {
325     object_count=2;
326     for (i_data=0; i_data<num_data;++i_data) {
327     if (! isEmpty(data_pp[i_data])) {
328     object_count++;
329     fprintf(fileHandle_p, "object \"%s\" class field\n",names_p[i_data]);
330     fprintf(fileHandle_p, "component \"positions\" value 1\n");
331     fprintf(fileHandle_p, "component \"connections\" value 2\n");
332     fprintf(fileHandle_p, "component \"data\" value %d\n",object_count);
333     }
334     }
335     }
336 jgs 82 fprintf(fileHandle_p, "end\n");
337     /* close the file */
338     fclose(fileHandle_p);
339 gross 1028 MEMFREE(isCellCentered);
340 jgs 82 return;
341     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26