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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26