/[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 1388 - (show annotations)
Fri Jan 11 07:45:58 2008 UTC (11 years, 9 months ago) by trankine
File MIME type: text/plain
File size: 11684 byte(s)
And get the *(&(*&(* name right
1
2 /* $Id$ */
3
4 /*******************************************************
5 *
6 * Copyright 2003-2007 by ACceSS MNRF
7 * Copyright 2007 by University of Queensland
8 *
9 * http://esscc.uq.edu.au
10 * Primary Business: Queensland, Australia
11 * Licensed under the Open Software License version 3.0
12 * http://www.opensource.org/licenses/osl-3.0.php
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 Finley_Mesh_saveDX(const char * filename_p, Finley_Mesh *mesh_p, const dim_t num_data,char* *names_p,escriptDataC* *data_pp) {
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 double* values,rtmp;
42 bool_t *isCellCentered=NULL;
43 Finley_ElementFile* elements=NULL;
44 ElementTypeId TypeId;
45 /* open the file and check handel */
46
47 /* if there is no mesh we just return */
48 if (mesh_p==NULL) return;
49 isCellCentered=MEMALLOC(num_data, bool_t);
50 if (Finley_checkPtr(isCellCentered)) return;
51
52 fileHandle_p = fopen(filename_p, "w");
53 if (fileHandle_p==NULL) {
54 sprintf(error_msg,"File %s could not be opened for writing.",filename_p);
55 MEMFREE(isCellCentered);
56 fclose(fileHandle_p);
57 Finley_setError(IO_ERROR,error_msg);
58 return;
59 }
60 /* find the mesh type to be written */
61 elementtype=FINLEY_UNKNOWN;
62 for (i_data=0;i_data<num_data;++i_data) {
63 if (! isEmpty(data_pp[i_data])) {
64 switch(getFunctionSpaceType(data_pp[i_data])) {
65 case FINLEY_REDUCED_NODES:
66 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
67 elementtype=FINLEY_ELEMENTS;
68 } else {
69 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data into single file.");
70 MEMFREE(isCellCentered);
71 fclose(fileHandle_p);
72 return;
73 }
74 isCellCentered[i_data]=FALSE;
75 break;
76 case FINLEY_ELEMENTS:
77 case FINLEY_REDUCED_ELEMENTS:
78 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
79 elementtype=FINLEY_ELEMENTS;
80 } else {
81 Finley_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]=TRUE;
87 break;
88 case FINLEY_FACE_ELEMENTS:
89 case FINLEY_REDUCED_FACE_ELEMENTS:
90 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_FACE_ELEMENTS) {
91 elementtype=FINLEY_FACE_ELEMENTS;
92 } else {
93 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data into single file.");
94 MEMFREE(isCellCentered);
95 fclose(fileHandle_p);
96 return;
97 }
98 isCellCentered[i_data]=TRUE;
99 break;
100 case FINLEY_POINTS:
101 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_POINTS) {
102 elementtype=FINLEY_POINTS;
103 } else {
104 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data into single file.");
105 MEMFREE(isCellCentered);
106 fclose(fileHandle_p);
107 return;
108 }
109 isCellCentered[i_data]=TRUE;
110 break;
111 case FINLEY_CONTACT_ELEMENTS_1:
112 case FINLEY_REDUCED_CONTACT_ELEMENTS_1:
113 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {
114 elementtype=FINLEY_CONTACT_ELEMENTS_1;
115 } else {
116 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data into single file.");
117 MEMFREE(isCellCentered);
118 fclose(fileHandle_p);
119 return;
120 }
121 isCellCentered[i_data]=TRUE;
122 break;
123 case FINLEY_CONTACT_ELEMENTS_2:
124 case FINLEY_REDUCED_CONTACT_ELEMENTS_2:
125 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {
126 elementtype=FINLEY_CONTACT_ELEMENTS_1;
127 } else {
128 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data into single file.");
129 MEMFREE(isCellCentered);
130 fclose(fileHandle_p);
131 return;
132 }
133 isCellCentered[i_data]=TRUE;
134 break;
135 default:
136 sprintf(error_msg,"saveDX: I don't know how to handle function space type %d",getFunctionSpaceType(data_pp[i_data]));
137 Finley_setError(TYPE_ERROR,error_msg);
138 MEMFREE(isCellCentered);
139 fclose(fileHandle_p);
140 return;
141 }
142 }
143 }
144 /* select number of points and the mesh component */
145 numPoints = mesh_p->Nodes->reducedNodesMapping->numTargets;
146 nDim = mesh_p->Nodes->numDim;
147 if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;
148 switch(elementtype) {
149 case FINLEY_ELEMENTS:
150 elements=mesh_p->Elements;
151 break;
152 case FINLEY_FACE_ELEMENTS:
153 elements=mesh_p->FaceElements;
154 break;
155 case FINLEY_POINTS:
156 elements=mesh_p->Points;
157 break;
158 case FINLEY_CONTACT_ELEMENTS_2:
159 elements=mesh_p->ContactElements;
160 break;
161 case FINLEY_CONTACT_ELEMENTS_1:
162 elements=mesh_p->ContactElements;
163 break;
164 }
165 if (elements==NULL) {
166 Finley_setError(SYSTEM_ERROR,"saveDX: undefined element file");
167 MEMFREE(isCellCentered);
168 fclose(fileHandle_p);
169 return;
170 }
171
172 /* map finley element type to DX element type */
173 TypeId = elements->ReferenceElement->Type->TypeId;
174 numDXNodesPerElement=0;
175 numCells = elements->numElements;
176 if (TypeId==Line2 || TypeId==Line3 || TypeId==Line4 ) {
177 numDXNodesPerElement=2;
178 resortIndex=resort[0];
179 strcpy(elemTypeStr, "lines");
180 } else if (TypeId==Tri3 || TypeId==Tri6 || TypeId==Tri9 || TypeId==Tri10 ) {
181 numDXNodesPerElement = 3;
182 resortIndex=resort[1];
183 strcpy(elemTypeStr, "triangles");
184 } else if (TypeId==Rec4 || TypeId==Rec8 || TypeId==Rec9 || TypeId==Rec12 || TypeId==Rec16 ) {
185 numDXNodesPerElement = 4;
186 resortIndex=resort[3];
187 strcpy(elemTypeStr, "quads");
188 } else if (TypeId==Tet4 || TypeId==Tet10 || TypeId==Tet16 ) {
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 } else {
197 sprintf(error_msg,"saveDX: Element type %s is not supported by DX",elements->ReferenceElement->Type->Name);
198 Finley_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 p=mesh_p->Nodes->reducedNodesMapping->map[i];
208 for (j = 0; j < nDim; j++) fprintf(fileHandle_p, " %g",mesh_p->Nodes->Coordinates[INDEX2(j, p, nDim)]);
209 fprintf(fileHandle_p, "\n");
210 }
211 /* connection table */
212 NN=elements->numNodes;
213 fprintf(fileHandle_p, "object 2 class array type int rank 1 shape %d items %d data follows\n",numDXNodesPerElement, numCells);
214 for (i = 0; i < numCells; i++) {
215 for (j = 0; j < numDXNodesPerElement; j++)
216 fprintf(fileHandle_p," %d",mesh_p->Nodes->reducedNodesMapping->target[elements->Nodes[INDEX2(resortIndex[j], i, NN)]]);
217 fprintf(fileHandle_p, "\n");
218 }
219 fprintf(fileHandle_p, "attribute \"element type\" string \"%s\"\n",elemTypeStr);
220 fprintf(fileHandle_p, "attribute \"ref\" string \"positions\"\n");
221
222 /* data */
223 object_count=2;
224 for (i_data =0 ;i_data<num_data;++i_data) {
225 if (! isEmpty(data_pp[i_data])) {
226 object_count++;
227 rank=getDataPointRank(data_pp[i_data]);
228 nComp=getDataPointSize(data_pp[i_data]);
229 fprintf(fileHandle_p, "object %d class array type float rank %d ",object_count,rank);
230 if (0 < rank) {
231 fprintf(fileHandle_p, "shape ");
232 for (i = 0; i < rank; i++) fprintf(fileHandle_p, "%d ", getDataPointShape(data_pp[i_data],i));
233 }
234 if (isCellCentered[i_data]) {
235 if (Finley_Assemble_reducedIntegrationOrder(data_pp[i_data])) {
236 numPointsPerSample=elements->ReferenceElementReducedOrder->numQuadNodes;
237 } else {
238 numPointsPerSample=elements->ReferenceElement->numQuadNodes;
239 }
240 if (numPointsPerSample>0) {
241 fprintf(fileHandle_p, "items %d data follows\n", numCells);
242 for (i=0;i<elements->numElements;i++) {
243 values=getSampleData(data_pp[i_data],i);
244 for (k=0;k<nComp;k++) {
245 if ( isExpanded(data_pp[i_data]) ) {
246 rtmp=0.;
247 for (j=0;j<numPointsPerSample;j++) rtmp+=values[INDEX2(k,j,nComp)];
248 fprintf(fileHandle_p, " %g", rtmp/numPointsPerSample);
249 } else {
250 fprintf(fileHandle_p, " %g", values[k]);
251 }
252 }
253 fprintf(fileHandle_p, "\n");
254 }
255 fprintf(fileHandle_p, "attribute \"dep\" string \"connections\"\n");
256 }
257 } else {
258 fprintf(fileHandle_p, "items %d data follows\n", numPoints);
259 for (i=0;i<numPoints;i++) {
260 values=getSampleData(data_pp[i_data],i);
261 for (k=0;k<nComp;k++) fprintf(fileHandle_p, " %g", values[k]);
262 fprintf(fileHandle_p, "\n");
263 }
264 fprintf(fileHandle_p, "attribute \"dep\" string \"positions\"\n");
265 }
266 }
267 }
268
269 /* and finish it up */
270 if (num_data==0) {
271 fprintf(fileHandle_p, "object %d class field\n",object_count+1);
272 fprintf(fileHandle_p, "component \"positions\" value 1\n");
273 fprintf(fileHandle_p, "component \"connections\" value 2\n");
274 } else {
275 object_count=2;
276 for (i_data=0; i_data<num_data;++i_data) {
277 if (! isEmpty(data_pp[i_data])) {
278 object_count++;
279 fprintf(fileHandle_p, "object \"%s\" class field\n",names_p[i_data]);
280 fprintf(fileHandle_p, "component \"positions\" value 1\n");
281 fprintf(fileHandle_p, "component \"connections\" value 2\n");
282 fprintf(fileHandle_p, "component \"data\" value %d\n",object_count);
283 }
284 }
285 }
286 fprintf(fileHandle_p, "end\n");
287 /* close the file */
288 fclose(fileHandle_p);
289 MEMFREE(isCellCentered);
290 return;
291 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26