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

Annotation of /trunk/dudley/src/Mesh_saveDX.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3221 - (hide annotations)
Wed Sep 29 01:00:21 2010 UTC (9 years, 2 months ago) by jfenwick
Original Path: branches/domexper/dudley/src/Mesh_saveDX.c
File MIME type: text/plain
File size: 10236 byte(s)
Comment stripping

1 jgs 82
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4 jfenwick 2881 * Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1811 * 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 jgs 82
14 ksteube 1811
15 jgs 82 /**************************************************************/
16    
17 ksteube 1312 /* writes data and mesh in an opendx file */
18     /* the input data needs to be cell centered or on reducedNodes */
19 jgs 82
20     /**************************************************************/
21    
22     #include "Mesh.h"
23 gross 1062 #include "Assemble.h"
24 jgs 82
25 jgs 150 /**************************************************************/
26    
27 jfenwick 3086 void Dudley_Mesh_saveDX(const char * filename_p, Dudley_Mesh *mesh_p, const dim_t num_data,char* *names_p,escriptDataC* *data_pp) {
28 gross 1028 char error_msg[LenErrorMsg_MAX], elemTypeStr[32];
29 jgs 82 /* 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 gross 1028 FILE * fileHandle_p = NULL;
38 ksteube 1312 int i,j,k,i_data, elementtype, numPoints = 0, nDim, *resortIndex=NULL, p,
39 gross 1062 numDXNodesPerElement=0, numCells, NN, object_count, rank, nComp, numPointsPerSample;
40 jfenwick 2271 __const double* values;
41     double rtmp;
42 gross 1028 bool_t *isCellCentered=NULL;
43 jfenwick 3086 Dudley_ElementFile* elements=NULL;
44 gross 1028 ElementTypeId TypeId;
45 jgs 82 /* open the file and check handel */
46 jgs 153
47 gross 1028 /* if there is no mesh we just return */
48     if (mesh_p==NULL) return;
49     isCellCentered=MEMALLOC(num_data, bool_t);
50 jfenwick 3086 if (Dudley_checkPtr(isCellCentered)) return;
51 gross 1028
52     fileHandle_p = fopen(filename_p, "w");
53 jgs 82 if (fileHandle_p==NULL) {
54 jgs 153 sprintf(error_msg,"File %s could not be opened for writing.",filename_p);
55 gross 1028 MEMFREE(isCellCentered);
56     fclose(fileHandle_p);
57 jfenwick 3086 Dudley_setError(IO_ERROR,error_msg);
58 jgs 153 return;
59 jgs 82 }
60 jgs 153 /* find the mesh type to be written */
61 jfenwick 3086 elementtype=DUDLEY_UNKNOWN;
62 jgs 153 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 jfenwick 3086 case DUDLEY_REDUCED_NODES:
66     if (elementtype==DUDLEY_UNKNOWN || elementtype==DUDLEY_ELEMENTS) {
67     elementtype=DUDLEY_ELEMENTS;
68 jgs 153 } else {
69 jfenwick 3086 Dudley_setError(TYPE_ERROR,"saveDX: cannot write given data into single file.");
70 gross 1028 MEMFREE(isCellCentered);
71     fclose(fileHandle_p);
72 jgs 153 return;
73     }
74     isCellCentered[i_data]=FALSE;
75     break;
76 jfenwick 3086 case DUDLEY_ELEMENTS:
77     case DUDLEY_REDUCED_ELEMENTS:
78     if (elementtype==DUDLEY_UNKNOWN || elementtype==DUDLEY_ELEMENTS) {
79     elementtype=DUDLEY_ELEMENTS;
80 jgs 153 } else {
81 jfenwick 3086 Dudley_setError(TYPE_ERROR,"saveDX: cannot write given data into single file.");
82 gross 1028 MEMFREE(isCellCentered);
83     fclose(fileHandle_p);
84 jgs 153 return;
85     }
86     isCellCentered[i_data]=TRUE;
87     break;
88 jfenwick 3086 case DUDLEY_FACE_ELEMENTS:
89     case DUDLEY_REDUCED_FACE_ELEMENTS:
90     if (elementtype==DUDLEY_UNKNOWN || elementtype==DUDLEY_FACE_ELEMENTS) {
91     elementtype=DUDLEY_FACE_ELEMENTS;
92 jgs 153 } else {
93 jfenwick 3086 Dudley_setError(TYPE_ERROR,"saveDX: cannot write given data into single file.");
94 gross 1028 MEMFREE(isCellCentered);
95     fclose(fileHandle_p);
96 jgs 153 return;
97     }
98     isCellCentered[i_data]=TRUE;
99     break;
100 jfenwick 3086 case DUDLEY_POINTS:
101     if (elementtype==DUDLEY_UNKNOWN || elementtype==DUDLEY_POINTS) {
102     elementtype=DUDLEY_POINTS;
103 jgs 153 } else {
104 jfenwick 3086 Dudley_setError(TYPE_ERROR,"saveDX: cannot write given data into single file.");
105 gross 1028 MEMFREE(isCellCentered);
106     fclose(fileHandle_p);
107 jgs 153 return;
108     }
109     isCellCentered[i_data]=TRUE;
110     break;
111     default:
112 gross 1361 sprintf(error_msg,"saveDX: I don't know how to handle function space type %d",getFunctionSpaceType(data_pp[i_data]));
113 jfenwick 3086 Dudley_setError(TYPE_ERROR,error_msg);
114 gross 1028 MEMFREE(isCellCentered);
115     fclose(fileHandle_p);
116 jgs 153 return;
117     }
118     }
119     }
120     /* select number of points and the mesh component */
121 ksteube 1312 numPoints = mesh_p->Nodes->reducedNodesMapping->numTargets;
122 gross 1028 nDim = mesh_p->Nodes->numDim;
123 jfenwick 3086 if (elementtype==DUDLEY_UNKNOWN) elementtype=DUDLEY_ELEMENTS;
124 jgs 153 switch(elementtype) {
125 jfenwick 3086 case DUDLEY_ELEMENTS:
126 jgs 153 elements=mesh_p->Elements;
127     break;
128 jfenwick 3086 case DUDLEY_FACE_ELEMENTS:
129 jgs 153 elements=mesh_p->FaceElements;
130     break;
131 jfenwick 3086 case DUDLEY_POINTS:
132 jgs 153 elements=mesh_p->Points;
133 jfenwick 3090 break; }
134 jgs 153 if (elements==NULL) {
135 jfenwick 3086 Dudley_setError(SYSTEM_ERROR,"saveDX: undefined element file");
136 gross 1028 MEMFREE(isCellCentered);
137     fclose(fileHandle_p);
138 jgs 153 return;
139     }
140    
141 jfenwick 3086 /* map dudley element type to DX element type */
142 jfenwick 3216 TypeId = elements->etype;
143 gross 1028 numDXNodesPerElement=0;
144     numCells = elements->numElements;
145 jfenwick 3126 if (TypeId==Line2 ) {
146 jgs 153 numDXNodesPerElement=2;
147     resortIndex=resort[0];
148     strcpy(elemTypeStr, "lines");
149 jfenwick 3126 } else if (TypeId==Tri3 ) {
150 jgs 153 numDXNodesPerElement = 3;
151     resortIndex=resort[1];
152     strcpy(elemTypeStr, "triangles");
153 jfenwick 3114 /* } else if (TypeId==Rec4 || TypeId==Rec8 || TypeId==Rec9 || TypeId==Rec12 || TypeId==Rec16 || TypeId==Rec9Macro ) {
154 jgs 153 numDXNodesPerElement = 4;
155     resortIndex=resort[3];
156 jfenwick 3114 strcpy(elemTypeStr, "quads");*/
157 jfenwick 3126 } else if (TypeId==Tet4 ) {
158 jgs 153 numDXNodesPerElement = 4;
159     resortIndex=resort[2];
160     strcpy(elemTypeStr, "tetrahedra");
161 jfenwick 3114 /* } else if (TypeId==Hex8 || TypeId==Hex20 || TypeId==Hex32 ) {
162 jgs 153 numDXNodesPerElement = 8;
163     resortIndex=resort[4];
164 jfenwick 3114 strcpy(elemTypeStr, "cubes");*/
165 jgs 153 } else {
166 jfenwick 3216 sprintf(error_msg,"saveDX: Element type %s is not supported by DX",elements->ename);
167 jfenwick 3086 Dudley_setError(VALUE_ERROR,error_msg);
168 gross 1028 MEMFREE(isCellCentered);
169     fclose(fileHandle_p);
170 jgs 153 return;
171     }
172    
173 jgs 82 /* positions */
174 ksteube 1312 fprintf(fileHandle_p, "object 1 class array type float rank 1 shape %d items %d data follows\n",nDim, numPoints);
175     for (i = 0; i < numPoints; i++) {
176     p=mesh_p->Nodes->reducedNodesMapping->map[i];
177     for (j = 0; j < nDim; j++) fprintf(fileHandle_p, " %g",mesh_p->Nodes->Coordinates[INDEX2(j, p, nDim)]);
178 jgs 82 fprintf(fileHandle_p, "\n");
179     }
180 jgs 153 /* connection table */
181 ksteube 1312 NN=elements->numNodes;
182 jgs 153 fprintf(fileHandle_p, "object 2 class array type int rank 1 shape %d items %d data follows\n",numDXNodesPerElement, numCells);
183     for (i = 0; i < numCells; i++) {
184 ksteube 1312 for (j = 0; j < numDXNodesPerElement; j++)
185     fprintf(fileHandle_p," %d",mesh_p->Nodes->reducedNodesMapping->target[elements->Nodes[INDEX2(resortIndex[j], i, NN)]]);
186 jgs 153 fprintf(fileHandle_p, "\n");
187     }
188     fprintf(fileHandle_p, "attribute \"element type\" string \"%s\"\n",elemTypeStr);
189     fprintf(fileHandle_p, "attribute \"ref\" string \"positions\"\n");
190 jgs 82
191     /* data */
192 gross 1028 object_count=2;
193 jgs 153 for (i_data =0 ;i_data<num_data;++i_data) {
194     if (! isEmpty(data_pp[i_data])) {
195     object_count++;
196 gross 1028 rank=getDataPointRank(data_pp[i_data]);
197     nComp=getDataPointSize(data_pp[i_data]);
198 jgs 153 fprintf(fileHandle_p, "object %d class array type float rank %d ",object_count,rank);
199     if (0 < rank) {
200     fprintf(fileHandle_p, "shape ");
201     for (i = 0; i < rank; i++) fprintf(fileHandle_p, "%d ", getDataPointShape(data_pp[i_data],i));
202     }
203     if (isCellCentered[i_data]) {
204 jfenwick 3086 if (Dudley_Assemble_reducedIntegrationOrder(data_pp[i_data])) {
205 jfenwick 3216 numPointsPerSample=(elements->numLocalDim==0)?0:1;
206 gross 1062 } else {
207 jfenwick 3216 numPointsPerSample=(elements->numLocalDim==0)?0:(elements->numLocalDim+1);
208 gross 1062 }
209 jgs 153 if (numPointsPerSample>0) {
210     fprintf(fileHandle_p, "items %d data follows\n", numCells);
211     for (i=0;i<elements->numElements;i++) {
212 jfenwick 2770 values=getSampleDataRO(data_pp[i_data],i);
213 jgs 153 for (k=0;k<nComp;k++) {
214 gross 903 if ( isExpanded(data_pp[i_data]) ) {
215     rtmp=0.;
216     for (j=0;j<numPointsPerSample;j++) rtmp+=values[INDEX2(k,j,nComp)];
217     fprintf(fileHandle_p, " %g", rtmp/numPointsPerSample);
218     } else {
219     fprintf(fileHandle_p, " %g", values[k]);
220     }
221 jgs 153 }
222     fprintf(fileHandle_p, "\n");
223     }
224     fprintf(fileHandle_p, "attribute \"dep\" string \"connections\"\n");
225     }
226     } else {
227 ksteube 1312 fprintf(fileHandle_p, "items %d data follows\n", numPoints);
228     for (i=0;i<numPoints;i++) {
229 jfenwick 2770 values=getSampleDataRO(data_pp[i_data],i);
230 ksteube 1312 for (k=0;k<nComp;k++) fprintf(fileHandle_p, " %g", values[k]);
231     fprintf(fileHandle_p, "\n");
232 jgs 82 }
233 jgs 153 fprintf(fileHandle_p, "attribute \"dep\" string \"positions\"\n");
234 jgs 82 }
235 jgs 153 }
236 jgs 82 }
237    
238     /* and finish it up */
239 jgs 153 if (num_data==0) {
240     fprintf(fileHandle_p, "object %d class field\n",object_count+1);
241     fprintf(fileHandle_p, "component \"positions\" value 1\n");
242     fprintf(fileHandle_p, "component \"connections\" value 2\n");
243     } else {
244     object_count=2;
245     for (i_data=0; i_data<num_data;++i_data) {
246     if (! isEmpty(data_pp[i_data])) {
247     object_count++;
248     fprintf(fileHandle_p, "object \"%s\" class field\n",names_p[i_data]);
249     fprintf(fileHandle_p, "component \"positions\" value 1\n");
250     fprintf(fileHandle_p, "component \"connections\" value 2\n");
251     fprintf(fileHandle_p, "component \"data\" value %d\n",object_count);
252     }
253     }
254     }
255 jgs 82 fprintf(fileHandle_p, "end\n");
256     /* close the file */
257     fclose(fileHandle_p);
258 gross 1028 MEMFREE(isCellCentered);
259 jgs 82 return;
260     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26