/[escript]/trunk/dudley/src/Mesh_saveDX.c
ViewVC logotype

Contents of /trunk/dudley/src/Mesh_saveDX.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3224 - (show annotations)
Wed Sep 29 05:19:37 2010 UTC (9 years, 1 month ago) by jfenwick
Original Path: branches/domexper/dudley/src/Mesh_saveDX.c
File MIME type: text/plain
File size: 9439 byte(s)
indent -linux -nce -i4 -bl -bli0 -l120

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 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
14 /**************************************************************/
15
16 /* writes data and mesh in an opendx file */
17 /* the input data needs to be cell centered or on reducedNodes */
18
19 /**************************************************************/
20
21 #include "Mesh.h"
22 #include "Assemble.h"
23
24 /**************************************************************/
25
26 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 ElementTypeId TypeId;
46 /* open the file and check handel */
47
48 /* 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
55 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 case DUDLEY_ELEMENTS:
148 elements = mesh_p->Elements;
149 break;
150 case DUDLEY_FACE_ELEMENTS:
151 elements = mesh_p->FaceElements;
152 break;
153 case DUDLEY_POINTS:
154 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
165 /* map dudley element type to DX element type */
166 TypeId = elements->etype;
167 numDXNodesPerElement = 0;
168 numCells = elements->numElements;
169 if (TypeId == Line2)
170 {
171 numDXNodesPerElement = 2;
172 resortIndex = resort[0];
173 strcpy(elemTypeStr, "lines");
174 }
175 else if (TypeId == Tri3)
176 {
177 numDXNodesPerElement = 3;
178 resortIndex = resort[1];
179 strcpy(elemTypeStr, "triangles");
180 /* } else if (TypeId==Rec4 || TypeId==Rec8 || TypeId==Rec9 || TypeId==Rec12 || TypeId==Rec16 || TypeId==Rec9Macro ) {
181 numDXNodesPerElement = 4;
182 resortIndex=resort[3];
183 strcpy(elemTypeStr, "quads");*/
184 }
185 else if (TypeId == Tet4)
186 {
187 numDXNodesPerElement = 4;
188 resortIndex = resort[2];
189 strcpy(elemTypeStr, "tetrahedra");
190 /* } else if (TypeId==Hex8 || TypeId==Hex20 || TypeId==Hex32 ) {
191 numDXNodesPerElement = 8;
192 resortIndex=resort[4];
193 strcpy(elemTypeStr, "cubes");*/
194 }
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
204 /* 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
227 /* 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
293 /* 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 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26