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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1062 - (show 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 /*
2 ************************************************************
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 */
12
13
14 /**************************************************************/
15
16 /* writes data and mesh in an opendx file */
17
18 /**************************************************************/
19
20 /* Author: Lutz Gross, gross@access.edu.au */
21 /* Version: $Id$ */
22
23 /**************************************************************/
24
25 #include "Mesh.h"
26 #include "Assemble.h"
27
28 /**************************************************************/
29
30 void Finley_Mesh_saveDX(const char * filename_p, Finley_Mesh *mesh_p, const dim_t num_data,char* *names_p,escriptDataC* *data_pp) {
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, nodetype, elementtype, numPoints = 0, nDim, *resortIndex=NULL,
42 numDXNodesPerElement=0, numCells, NN, object_count, rank, nComp, numPointsPerSample;
43 double* values,rtmp;
44 bool_t *isCellCentered=NULL;
45 Finley_ElementFile* elements=NULL;
46 ElementTypeId TypeId;
47 /* open the file and check handel */
48
49 /* 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 if (fileHandle_p==NULL) {
56 sprintf(error_msg,"File %s could not be opened for writing.",filename_p);
57 MEMFREE(isCellCentered);
58 fclose(fileHandle_p);
59 Finley_setError(IO_ERROR,error_msg);
60 return;
61 }
62 /* find the mesh type to be written */
63 nodetype=FINLEY_DEGREES_OF_FREEDOM;
64 elementtype=FINLEY_UNKNOWN;
65 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 MEMFREE(isCellCentered);
75 fclose(fileHandle_p);
76 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 MEMFREE(isCellCentered);
87 fclose(fileHandle_p);
88 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 MEMFREE(isCellCentered);
99 fclose(fileHandle_p);
100 return;
101 }
102 isCellCentered[i_data]=FALSE;
103 break;
104 case FINLEY_ELEMENTS:
105 case FINLEY_REDUCED_ELEMENTS:
106 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 MEMFREE(isCellCentered);
112 fclose(fileHandle_p);
113 return;
114 }
115 isCellCentered[i_data]=TRUE;
116 break;
117 case FINLEY_FACE_ELEMENTS:
118 case FINLEY_REDUCED_FACE_ELEMENTS:
119 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 MEMFREE(isCellCentered);
125 fclose(fileHandle_p);
126 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 MEMFREE(isCellCentered);
137 fclose(fileHandle_p);
138 return;
139 }
140 isCellCentered[i_data]=TRUE;
141 break;
142 case FINLEY_CONTACT_ELEMENTS_1:
143 case FINLEY_REDUCED_CONTACT_ELEMENTS_1:
144 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 MEMFREE(isCellCentered);
150 fclose(fileHandle_p);
151 return;
152 }
153 isCellCentered[i_data]=TRUE;
154 break;
155 case FINLEY_CONTACT_ELEMENTS_2:
156 case FINLEY_REDUCED_CONTACT_ELEMENTS_2:
157 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 MEMFREE(isCellCentered);
163 fclose(fileHandle_p);
164 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 MEMFREE(isCellCentered);
172 fclose(fileHandle_p);
173 return;
174 }
175 }
176 }
177 /* select number of points and the mesh component */
178 numPoints = mesh_p->Nodes->numNodes;
179 nDim = mesh_p->Nodes->numDim;
180 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 case FINLEY_CONTACT_ELEMENTS_2:
197 elements=mesh_p->ContactElements;
198 break;
199 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 MEMFREE(isCellCentered);
206 fclose(fileHandle_p);
207 return;
208 }
209
210 /* map finley element type to DX element type */
211 TypeId = elements->ReferenceElement->Type->TypeId;
212 numDXNodesPerElement=0;
213 numCells = elements->numElements;
214 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 MEMFREE(isCellCentered);
238 fclose(fileHandle_p);
239 return;
240 }
241
242 /* positions */
243 fprintf(fileHandle_p, "object 1 class array type float rank 1 shape %d items %d data follows\n",nDim, mesh_p->Nodes->reducedNumNodes);
244 for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
245 if (mesh_p->Nodes->toReduced[i]>=0) {
246 for (j = 0; j < nDim; j++) fprintf(fileHandle_p, " %g",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);
247 fprintf(fileHandle_p, "\n");
248 }
249 }
250 /* connection table */
251 NN=elements->ReferenceElement->Type->numNodes;
252 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
260 /* data */
261 object_count=2;
262 for (i_data =0 ;i_data<num_data;++i_data) {
263 if (! isEmpty(data_pp[i_data])) {
264 object_count++;
265 rank=getDataPointRank(data_pp[i_data]);
266 nComp=getDataPointSize(data_pp[i_data]);
267 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 if (Finley_Assemble_reducedIntegrationOrder(data_pp[i_data])) {
274 numPointsPerSample=elements->ReferenceElementReducedOrder->numQuadNodes;
275 } else {
276 numPointsPerSample=elements->ReferenceElement->numQuadNodes;
277 }
278 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 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 }
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 }
313 }
314 fprintf(fileHandle_p, "attribute \"dep\" string \"positions\"\n");
315 }
316 }
317 }
318
319 /* and finish it up */
320 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 fprintf(fileHandle_p, "end\n");
337 /* close the file */
338 fclose(fileHandle_p);
339 MEMFREE(isCellCentered);
340 return;
341 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26