/[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 3247 - (hide annotations)
Wed Oct 6 05:53:06 2010 UTC (9 years, 2 months ago) by caltinay
Original Path: branches/domexper/dudley/src/Mesh_saveDX.c
File MIME type: text/plain
File size: 9467 byte(s)
Fixed name clashes between dudley and finley so both can be used
simultaneously.

1 jgs 82
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4 jfenwick 2881 * Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1811 * 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     /**************************************************************/
15    
16 ksteube 1312 /* writes data and mesh in an opendx file */
17     /* the input data needs to be cell centered or on reducedNodes */
18 jgs 82
19     /**************************************************************/
20    
21     #include "Mesh.h"
22 gross 1062 #include "Assemble.h"
23 jgs 82
24 jgs 150 /**************************************************************/
25    
26 jfenwick 3224 void Dudley_Mesh_saveDX(const char *filename_p, Dudley_Mesh * mesh_p, const dim_t num_data, char **names_p,
27     escriptDataC * *data_pp)
28     {
29     char error_msg[LenErrorMsg_MAX], elemTypeStr[32];
30     /* 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     FILE *fileHandle_p = NULL;
39     int i, j, k, i_data, elementtype, numPoints = 0, nDim, *resortIndex = NULL, p,
40     numDXNodesPerElement = 0, numCells, NN, object_count, rank, nComp, numPointsPerSample;
41     __const double *values;
42     double rtmp;
43     bool_t *isCellCentered = NULL;
44     Dudley_ElementFile *elements = NULL;
45 caltinay 3247 Dudley_ElementTypeId TypeId;
46 jfenwick 3224 /* open the file and check handel */
47 jgs 153
48 jfenwick 3224 /* if there is no mesh we just return */
49     if (mesh_p == NULL)
50     return;
51     isCellCentered = MEMALLOC(num_data, bool_t);
52     if (Dudley_checkPtr(isCellCentered))
53     return;
54 gross 1028
55 jfenwick 3224 fileHandle_p = fopen(filename_p, "w");
56     if (fileHandle_p == NULL)
57     {
58     sprintf(error_msg, "File %s could not be opened for writing.", filename_p);
59     MEMFREE(isCellCentered);
60     fclose(fileHandle_p);
61     Dudley_setError(IO_ERROR, error_msg);
62     return;
63     }
64     /* find the mesh type to be written */
65     elementtype = DUDLEY_UNKNOWN;
66     for (i_data = 0; i_data < num_data; ++i_data)
67     {
68     if (!isEmpty(data_pp[i_data]))
69     {
70     switch (getFunctionSpaceType(data_pp[i_data]))
71     {
72     case DUDLEY_REDUCED_NODES:
73     if (elementtype == DUDLEY_UNKNOWN || elementtype == DUDLEY_ELEMENTS)
74     {
75     elementtype = DUDLEY_ELEMENTS;
76     }
77     else
78     {
79     Dudley_setError(TYPE_ERROR, "saveDX: cannot write given data into single file.");
80     MEMFREE(isCellCentered);
81     fclose(fileHandle_p);
82     return;
83     }
84     isCellCentered[i_data] = FALSE;
85     break;
86     case DUDLEY_ELEMENTS:
87     case DUDLEY_REDUCED_ELEMENTS:
88     if (elementtype == DUDLEY_UNKNOWN || elementtype == DUDLEY_ELEMENTS)
89     {
90     elementtype = DUDLEY_ELEMENTS;
91     }
92     else
93     {
94     Dudley_setError(TYPE_ERROR, "saveDX: cannot write given data into single file.");
95     MEMFREE(isCellCentered);
96     fclose(fileHandle_p);
97     return;
98     }
99     isCellCentered[i_data] = TRUE;
100     break;
101     case DUDLEY_FACE_ELEMENTS:
102     case DUDLEY_REDUCED_FACE_ELEMENTS:
103     if (elementtype == DUDLEY_UNKNOWN || elementtype == DUDLEY_FACE_ELEMENTS)
104     {
105     elementtype = DUDLEY_FACE_ELEMENTS;
106     }
107     else
108     {
109     Dudley_setError(TYPE_ERROR, "saveDX: cannot write given data into single file.");
110     MEMFREE(isCellCentered);
111     fclose(fileHandle_p);
112     return;
113     }
114     isCellCentered[i_data] = TRUE;
115     break;
116     case DUDLEY_POINTS:
117     if (elementtype == DUDLEY_UNKNOWN || elementtype == DUDLEY_POINTS)
118     {
119     elementtype = DUDLEY_POINTS;
120     }
121     else
122     {
123     Dudley_setError(TYPE_ERROR, "saveDX: cannot write given data into single file.");
124     MEMFREE(isCellCentered);
125     fclose(fileHandle_p);
126     return;
127     }
128     isCellCentered[i_data] = TRUE;
129     break;
130     default:
131     sprintf(error_msg, "saveDX: I don't know how to handle function space type %d",
132     getFunctionSpaceType(data_pp[i_data]));
133     Dudley_setError(TYPE_ERROR, error_msg);
134     MEMFREE(isCellCentered);
135     fclose(fileHandle_p);
136     return;
137     }
138     }
139     }
140     /* select number of points and the mesh component */
141     numPoints = mesh_p->Nodes->reducedNodesMapping->numTargets;
142     nDim = mesh_p->Nodes->numDim;
143     if (elementtype == DUDLEY_UNKNOWN)
144     elementtype = DUDLEY_ELEMENTS;
145     switch (elementtype)
146     {
147 jfenwick 3086 case DUDLEY_ELEMENTS:
148 jfenwick 3224 elements = mesh_p->Elements;
149     break;
150 jfenwick 3086 case DUDLEY_FACE_ELEMENTS:
151 jfenwick 3224 elements = mesh_p->FaceElements;
152     break;
153 jfenwick 3086 case DUDLEY_POINTS:
154 jfenwick 3224 elements = mesh_p->Points;
155     break;
156     }
157     if (elements == NULL)
158     {
159     Dudley_setError(SYSTEM_ERROR, "saveDX: undefined element file");
160     MEMFREE(isCellCentered);
161     fclose(fileHandle_p);
162     return;
163     }
164 jgs 153
165 jfenwick 3224 /* map dudley element type to DX element type */
166     TypeId = elements->etype;
167     numDXNodesPerElement = 0;
168     numCells = elements->numElements;
169 caltinay 3247 if (TypeId == Dudley_Line2)
170 jfenwick 3224 {
171     numDXNodesPerElement = 2;
172     resortIndex = resort[0];
173     strcpy(elemTypeStr, "lines");
174     }
175 caltinay 3247 else if (TypeId == Dudley_Tri3)
176 jfenwick 3224 {
177     numDXNodesPerElement = 3;
178     resortIndex = resort[1];
179     strcpy(elemTypeStr, "triangles");
180 jfenwick 3114 /* } else if (TypeId==Rec4 || TypeId==Rec8 || TypeId==Rec9 || TypeId==Rec12 || TypeId==Rec16 || TypeId==Rec9Macro ) {
181 jgs 153 numDXNodesPerElement = 4;
182     resortIndex=resort[3];
183 jfenwick 3114 strcpy(elemTypeStr, "quads");*/
184 jfenwick 3224 }
185 caltinay 3247 else if (TypeId == Dudley_Tet4)
186 jfenwick 3224 {
187     numDXNodesPerElement = 4;
188     resortIndex = resort[2];
189     strcpy(elemTypeStr, "tetrahedra");
190 jfenwick 3114 /* } else if (TypeId==Hex8 || TypeId==Hex20 || TypeId==Hex32 ) {
191 jgs 153 numDXNodesPerElement = 8;
192     resortIndex=resort[4];
193 jfenwick 3114 strcpy(elemTypeStr, "cubes");*/
194 jfenwick 3224 }
195     else
196     {
197     sprintf(error_msg, "saveDX: Element type %s is not supported by DX", elements->ename);
198     Dudley_setError(VALUE_ERROR, error_msg);
199     MEMFREE(isCellCentered);
200     fclose(fileHandle_p);
201     return;
202     }
203 jgs 153
204 jfenwick 3224 /* positions */
205     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     {
208     p = mesh_p->Nodes->reducedNodesMapping->map[i];
209     for (j = 0; j < nDim; j++)
210     fprintf(fileHandle_p, " %g", mesh_p->Nodes->Coordinates[INDEX2(j, p, nDim)]);
211     fprintf(fileHandle_p, "\n");
212     }
213     /* connection table */
214     NN = elements->numNodes;
215     fprintf(fileHandle_p, "object 2 class array type int rank 1 shape %d items %d data follows\n", numDXNodesPerElement,
216     numCells);
217     for (i = 0; i < numCells; i++)
218     {
219     for (j = 0; j < numDXNodesPerElement; j++)
220     fprintf(fileHandle_p, " %d",
221     mesh_p->Nodes->reducedNodesMapping->target[elements->Nodes[INDEX2(resortIndex[j], i, NN)]]);
222     fprintf(fileHandle_p, "\n");
223     }
224     fprintf(fileHandle_p, "attribute \"element type\" string \"%s\"\n", elemTypeStr);
225     fprintf(fileHandle_p, "attribute \"ref\" string \"positions\"\n");
226 jgs 82
227 jfenwick 3224 /* data */
228     object_count = 2;
229     for (i_data = 0; i_data < num_data; ++i_data)
230     {
231     if (!isEmpty(data_pp[i_data]))
232     {
233     object_count++;
234     rank = getDataPointRank(data_pp[i_data]);
235     nComp = getDataPointSize(data_pp[i_data]);
236     fprintf(fileHandle_p, "object %d class array type float rank %d ", object_count, rank);
237     if (0 < rank)
238     {
239     fprintf(fileHandle_p, "shape ");
240     for (i = 0; i < rank; i++)
241     fprintf(fileHandle_p, "%d ", getDataPointShape(data_pp[i_data], i));
242     }
243     if (isCellCentered[i_data])
244     {
245     if (Dudley_Assemble_reducedIntegrationOrder(data_pp[i_data]))
246     {
247     numPointsPerSample = (elements->numLocalDim == 0) ? 0 : 1;
248     }
249     else
250     {
251     numPointsPerSample = (elements->numLocalDim == 0) ? 0 : (elements->numLocalDim + 1);
252     }
253     if (numPointsPerSample > 0)
254     {
255     fprintf(fileHandle_p, "items %d data follows\n", numCells);
256     for (i = 0; i < elements->numElements; i++)
257     {
258     values = getSampleDataRO(data_pp[i_data], i);
259     for (k = 0; k < nComp; k++)
260     {
261     if (isExpanded(data_pp[i_data]))
262     {
263     rtmp = 0.;
264     for (j = 0; j < numPointsPerSample; j++)
265     rtmp += values[INDEX2(k, j, nComp)];
266     fprintf(fileHandle_p, " %g", rtmp / numPointsPerSample);
267     }
268     else
269     {
270     fprintf(fileHandle_p, " %g", values[k]);
271     }
272     }
273     fprintf(fileHandle_p, "\n");
274     }
275     fprintf(fileHandle_p, "attribute \"dep\" string \"connections\"\n");
276     }
277     }
278     else
279     {
280     fprintf(fileHandle_p, "items %d data follows\n", numPoints);
281     for (i = 0; i < numPoints; i++)
282     {
283     values = getSampleDataRO(data_pp[i_data], i);
284     for (k = 0; k < nComp; k++)
285     fprintf(fileHandle_p, " %g", values[k]);
286     fprintf(fileHandle_p, "\n");
287     }
288     fprintf(fileHandle_p, "attribute \"dep\" string \"positions\"\n");
289     }
290     }
291     }
292 jgs 82
293 jfenwick 3224 /* and finish it up */
294     if (num_data == 0)
295     {
296     fprintf(fileHandle_p, "object %d class field\n", object_count + 1);
297     fprintf(fileHandle_p, "component \"positions\" value 1\n");
298     fprintf(fileHandle_p, "component \"connections\" value 2\n");
299     }
300     else
301     {
302     object_count = 2;
303     for (i_data = 0; i_data < num_data; ++i_data)
304     {
305     if (!isEmpty(data_pp[i_data]))
306     {
307     object_count++;
308     fprintf(fileHandle_p, "object \"%s\" class field\n", names_p[i_data]);
309     fprintf(fileHandle_p, "component \"positions\" value 1\n");
310     fprintf(fileHandle_p, "component \"connections\" value 2\n");
311     fprintf(fileHandle_p, "component \"data\" value %d\n", object_count);
312     }
313     }
314     }
315     fprintf(fileHandle_p, "end\n");
316     /* close the file */
317     fclose(fileHandle_p);
318     MEMFREE(isCellCentered);
319     return;
320 jgs 82 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26