/[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 153 - (show annotations)
Tue Oct 25 01:51:20 2005 UTC (14 years, 1 month ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Mesh_saveDX.c
File MIME type: text/plain
File size: 13202 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-10-25

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26