/[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 1811 - (show annotations)
Thu Sep 25 23:11:13 2008 UTC (10 years, 11 months ago) by ksteube
File MIME type: text/plain
File size: 11649 byte(s)
Copyright updated in all files

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2008 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
17 /* writes data and mesh in an opendx file */
18 /* the input data needs to be cell centered or on reducedNodes */
19
20 /**************************************************************/
21
22 #include "Mesh.h"
23 #include "Assemble.h"
24
25 /**************************************************************/
26
27 void Finley_Mesh_saveDX(const char * filename_p, Finley_Mesh *mesh_p, const dim_t num_data,char* *names_p,escriptDataC* *data_pp) {
28 char error_msg[LenErrorMsg_MAX], elemTypeStr[32];
29 /* some tables needed for reordering */
30 int resort[6][9]={
31 {0,1}, /* line */
32 {0,1,2}, /* triagle */
33 {0,1,2,3}, /* tetrahedron */
34 {0,3,1,2}, /* quadrilateral */
35 {3,0,7,4,2,1,6,5}, /* hexahedron */
36 };
37 FILE * fileHandle_p = NULL;
38 int i,j,k,i_data, elementtype, numPoints = 0, nDim, *resortIndex=NULL, p,
39 numDXNodesPerElement=0, numCells, NN, object_count, rank, nComp, numPointsPerSample;
40 double* values,rtmp;
41 bool_t *isCellCentered=NULL;
42 Finley_ElementFile* elements=NULL;
43 ElementTypeId TypeId;
44 /* open the file and check handel */
45
46 /* if there is no mesh we just return */
47 if (mesh_p==NULL) return;
48 isCellCentered=MEMALLOC(num_data, bool_t);
49 if (Finley_checkPtr(isCellCentered)) return;
50
51 fileHandle_p = fopen(filename_p, "w");
52 if (fileHandle_p==NULL) {
53 sprintf(error_msg,"File %s could not be opened for writing.",filename_p);
54 MEMFREE(isCellCentered);
55 fclose(fileHandle_p);
56 Finley_setError(IO_ERROR,error_msg);
57 return;
58 }
59 /* find the mesh type to be written */
60 elementtype=FINLEY_UNKNOWN;
61 for (i_data=0;i_data<num_data;++i_data) {
62 if (! isEmpty(data_pp[i_data])) {
63 switch(getFunctionSpaceType(data_pp[i_data])) {
64 case FINLEY_REDUCED_NODES:
65 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
66 elementtype=FINLEY_ELEMENTS;
67 } else {
68 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data into single file.");
69 MEMFREE(isCellCentered);
70 fclose(fileHandle_p);
71 return;
72 }
73 isCellCentered[i_data]=FALSE;
74 break;
75 case FINLEY_ELEMENTS:
76 case FINLEY_REDUCED_ELEMENTS:
77 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
78 elementtype=FINLEY_ELEMENTS;
79 } else {
80 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data into single file.");
81 MEMFREE(isCellCentered);
82 fclose(fileHandle_p);
83 return;
84 }
85 isCellCentered[i_data]=TRUE;
86 break;
87 case FINLEY_FACE_ELEMENTS:
88 case FINLEY_REDUCED_FACE_ELEMENTS:
89 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_FACE_ELEMENTS) {
90 elementtype=FINLEY_FACE_ELEMENTS;
91 } else {
92 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data into single file.");
93 MEMFREE(isCellCentered);
94 fclose(fileHandle_p);
95 return;
96 }
97 isCellCentered[i_data]=TRUE;
98 break;
99 case FINLEY_POINTS:
100 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_POINTS) {
101 elementtype=FINLEY_POINTS;
102 } else {
103 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data into single file.");
104 MEMFREE(isCellCentered);
105 fclose(fileHandle_p);
106 return;
107 }
108 isCellCentered[i_data]=TRUE;
109 break;
110 case FINLEY_CONTACT_ELEMENTS_1:
111 case FINLEY_REDUCED_CONTACT_ELEMENTS_1:
112 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {
113 elementtype=FINLEY_CONTACT_ELEMENTS_1;
114 } else {
115 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data into single file.");
116 MEMFREE(isCellCentered);
117 fclose(fileHandle_p);
118 return;
119 }
120 isCellCentered[i_data]=TRUE;
121 break;
122 case FINLEY_CONTACT_ELEMENTS_2:
123 case FINLEY_REDUCED_CONTACT_ELEMENTS_2:
124 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {
125 elementtype=FINLEY_CONTACT_ELEMENTS_1;
126 } else {
127 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data into single file.");
128 MEMFREE(isCellCentered);
129 fclose(fileHandle_p);
130 return;
131 }
132 isCellCentered[i_data]=TRUE;
133 break;
134 default:
135 sprintf(error_msg,"saveDX: I don't know how to handle function space type %d",getFunctionSpaceType(data_pp[i_data]));
136 Finley_setError(TYPE_ERROR,error_msg);
137 MEMFREE(isCellCentered);
138 fclose(fileHandle_p);
139 return;
140 }
141 }
142 }
143 /* select number of points and the mesh component */
144 numPoints = mesh_p->Nodes->reducedNodesMapping->numTargets;
145 nDim = mesh_p->Nodes->numDim;
146 if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;
147 switch(elementtype) {
148 case FINLEY_ELEMENTS:
149 elements=mesh_p->Elements;
150 break;
151 case FINLEY_FACE_ELEMENTS:
152 elements=mesh_p->FaceElements;
153 break;
154 case FINLEY_POINTS:
155 elements=mesh_p->Points;
156 break;
157 case FINLEY_CONTACT_ELEMENTS_2:
158 elements=mesh_p->ContactElements;
159 break;
160 case FINLEY_CONTACT_ELEMENTS_1:
161 elements=mesh_p->ContactElements;
162 break;
163 }
164 if (elements==NULL) {
165 Finley_setError(SYSTEM_ERROR,"saveDX: undefined element file");
166 MEMFREE(isCellCentered);
167 fclose(fileHandle_p);
168 return;
169 }
170
171 /* map finley element type to DX element type */
172 TypeId = elements->ReferenceElement->Type->TypeId;
173 numDXNodesPerElement=0;
174 numCells = elements->numElements;
175 if (TypeId==Line2 || TypeId==Line3 || TypeId==Line4 ) {
176 numDXNodesPerElement=2;
177 resortIndex=resort[0];
178 strcpy(elemTypeStr, "lines");
179 } else if (TypeId==Tri3 || TypeId==Tri6 || TypeId==Tri9 || TypeId==Tri10 ) {
180 numDXNodesPerElement = 3;
181 resortIndex=resort[1];
182 strcpy(elemTypeStr, "triangles");
183 } else if (TypeId==Rec4 || TypeId==Rec8 || TypeId==Rec9 || TypeId==Rec12 || TypeId==Rec16 ) {
184 numDXNodesPerElement = 4;
185 resortIndex=resort[3];
186 strcpy(elemTypeStr, "quads");
187 } else if (TypeId==Tet4 || TypeId==Tet10 || TypeId==Tet16 ) {
188 numDXNodesPerElement = 4;
189 resortIndex=resort[2];
190 strcpy(elemTypeStr, "tetrahedra");
191 } else if (TypeId==Hex8 || TypeId==Hex20 || TypeId==Hex32 ) {
192 numDXNodesPerElement = 8;
193 resortIndex=resort[4];
194 strcpy(elemTypeStr, "cubes");
195 } else {
196 sprintf(error_msg,"saveDX: Element type %s is not supported by DX",elements->ReferenceElement->Type->Name);
197 Finley_setError(VALUE_ERROR,error_msg);
198 MEMFREE(isCellCentered);
199 fclose(fileHandle_p);
200 return;
201 }
202
203 /* positions */
204 fprintf(fileHandle_p, "object 1 class array type float rank 1 shape %d items %d data follows\n",nDim, numPoints);
205 for (i = 0; i < numPoints; i++) {
206 p=mesh_p->Nodes->reducedNodesMapping->map[i];
207 for (j = 0; j < nDim; j++) fprintf(fileHandle_p, " %g",mesh_p->Nodes->Coordinates[INDEX2(j, p, nDim)]);
208 fprintf(fileHandle_p, "\n");
209 }
210 /* connection table */
211 NN=elements->numNodes;
212 fprintf(fileHandle_p, "object 2 class array type int rank 1 shape %d items %d data follows\n",numDXNodesPerElement, numCells);
213 for (i = 0; i < numCells; i++) {
214 for (j = 0; j < numDXNodesPerElement; j++)
215 fprintf(fileHandle_p," %d",mesh_p->Nodes->reducedNodesMapping->target[elements->Nodes[INDEX2(resortIndex[j], i, NN)]]);
216 fprintf(fileHandle_p, "\n");
217 }
218 fprintf(fileHandle_p, "attribute \"element type\" string \"%s\"\n",elemTypeStr);
219 fprintf(fileHandle_p, "attribute \"ref\" string \"positions\"\n");
220
221 /* data */
222 object_count=2;
223 for (i_data =0 ;i_data<num_data;++i_data) {
224 if (! isEmpty(data_pp[i_data])) {
225 object_count++;
226 rank=getDataPointRank(data_pp[i_data]);
227 nComp=getDataPointSize(data_pp[i_data]);
228 fprintf(fileHandle_p, "object %d class array type float rank %d ",object_count,rank);
229 if (0 < rank) {
230 fprintf(fileHandle_p, "shape ");
231 for (i = 0; i < rank; i++) fprintf(fileHandle_p, "%d ", getDataPointShape(data_pp[i_data],i));
232 }
233 if (isCellCentered[i_data]) {
234 if (Finley_Assemble_reducedIntegrationOrder(data_pp[i_data])) {
235 numPointsPerSample=elements->ReferenceElementReducedOrder->numQuadNodes;
236 } else {
237 numPointsPerSample=elements->ReferenceElement->numQuadNodes;
238 }
239 if (numPointsPerSample>0) {
240 fprintf(fileHandle_p, "items %d data follows\n", numCells);
241 for (i=0;i<elements->numElements;i++) {
242 values=getSampleData(data_pp[i_data],i);
243 for (k=0;k<nComp;k++) {
244 if ( isExpanded(data_pp[i_data]) ) {
245 rtmp=0.;
246 for (j=0;j<numPointsPerSample;j++) rtmp+=values[INDEX2(k,j,nComp)];
247 fprintf(fileHandle_p, " %g", rtmp/numPointsPerSample);
248 } else {
249 fprintf(fileHandle_p, " %g", values[k]);
250 }
251 }
252 fprintf(fileHandle_p, "\n");
253 }
254 fprintf(fileHandle_p, "attribute \"dep\" string \"connections\"\n");
255 }
256 } else {
257 fprintf(fileHandle_p, "items %d data follows\n", numPoints);
258 for (i=0;i<numPoints;i++) {
259 values=getSampleData(data_pp[i_data],i);
260 for (k=0;k<nComp;k++) fprintf(fileHandle_p, " %g", values[k]);
261 fprintf(fileHandle_p, "\n");
262 }
263 fprintf(fileHandle_p, "attribute \"dep\" string \"positions\"\n");
264 }
265 }
266 }
267
268 /* and finish it up */
269 if (num_data==0) {
270 fprintf(fileHandle_p, "object %d class field\n",object_count+1);
271 fprintf(fileHandle_p, "component \"positions\" value 1\n");
272 fprintf(fileHandle_p, "component \"connections\" value 2\n");
273 } else {
274 object_count=2;
275 for (i_data=0; i_data<num_data;++i_data) {
276 if (! isEmpty(data_pp[i_data])) {
277 object_count++;
278 fprintf(fileHandle_p, "object \"%s\" class field\n",names_p[i_data]);
279 fprintf(fileHandle_p, "component \"positions\" value 1\n");
280 fprintf(fileHandle_p, "component \"connections\" value 2\n");
281 fprintf(fileHandle_p, "component \"data\" value %d\n",object_count);
282 }
283 }
284 }
285 fprintf(fileHandle_p, "end\n");
286 /* close the file */
287 fclose(fileHandle_p);
288 MEMFREE(isCellCentered);
289 return;
290 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26