/[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 3086 - (hide annotations)
Thu Aug 5 05:07:58 2010 UTC (9 years, 3 months ago) by jfenwick
Original Path: branches/domexper/dudley/src/Mesh_saveDX.c
File MIME type: text/plain
File size: 11859 byte(s)
Another pass at removing finley

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 jfenwick 3086 case DUDLEY_CONTACT_ELEMENTS_1:
112     case DUDLEY_REDUCED_CONTACT_ELEMENTS_1:
113     if (elementtype==DUDLEY_UNKNOWN || elementtype==DUDLEY_CONTACT_ELEMENTS_1) {
114     elementtype=DUDLEY_CONTACT_ELEMENTS_1;
115 jgs 153 } else {
116 jfenwick 3086 Dudley_setError(TYPE_ERROR,"saveDX: cannot write given data into single file.");
117 gross 1028 MEMFREE(isCellCentered);
118     fclose(fileHandle_p);
119 jgs 153 return;
120     }
121     isCellCentered[i_data]=TRUE;
122     break;
123 jfenwick 3086 case DUDLEY_CONTACT_ELEMENTS_2:
124     case DUDLEY_REDUCED_CONTACT_ELEMENTS_2:
125     if (elementtype==DUDLEY_UNKNOWN || elementtype==DUDLEY_CONTACT_ELEMENTS_1) {
126     elementtype=DUDLEY_CONTACT_ELEMENTS_1;
127 jgs 153 } else {
128 jfenwick 3086 Dudley_setError(TYPE_ERROR,"saveDX: cannot write given data into single file.");
129 gross 1028 MEMFREE(isCellCentered);
130     fclose(fileHandle_p);
131 jgs 153 return;
132     }
133     isCellCentered[i_data]=TRUE;
134     break;
135     default:
136 gross 1361 sprintf(error_msg,"saveDX: I don't know how to handle function space type %d",getFunctionSpaceType(data_pp[i_data]));
137 jfenwick 3086 Dudley_setError(TYPE_ERROR,error_msg);
138 gross 1028 MEMFREE(isCellCentered);
139     fclose(fileHandle_p);
140 jgs 153 return;
141     }
142     }
143     }
144     /* select number of points and the mesh component */
145 ksteube 1312 numPoints = mesh_p->Nodes->reducedNodesMapping->numTargets;
146 gross 1028 nDim = mesh_p->Nodes->numDim;
147 jfenwick 3086 if (elementtype==DUDLEY_UNKNOWN) elementtype=DUDLEY_ELEMENTS;
148 jgs 153 switch(elementtype) {
149 jfenwick 3086 case DUDLEY_ELEMENTS:
150 jgs 153 elements=mesh_p->Elements;
151     break;
152 jfenwick 3086 case DUDLEY_FACE_ELEMENTS:
153 jgs 153 elements=mesh_p->FaceElements;
154     break;
155 jfenwick 3086 case DUDLEY_POINTS:
156 jgs 153 elements=mesh_p->Points;
157     break;
158 jfenwick 3086 case DUDLEY_CONTACT_ELEMENTS_2:
159 gross 1062 elements=mesh_p->ContactElements;
160     break;
161 jfenwick 3086 case DUDLEY_CONTACT_ELEMENTS_1:
162 jgs 153 elements=mesh_p->ContactElements;
163     break;
164     }
165     if (elements==NULL) {
166 jfenwick 3086 Dudley_setError(SYSTEM_ERROR,"saveDX: undefined element file");
167 gross 1028 MEMFREE(isCellCentered);
168     fclose(fileHandle_p);
169 jgs 153 return;
170     }
171    
172 jfenwick 3086 /* map dudley element type to DX element type */
173 gross 2748 TypeId = elements->referenceElementSet->referenceElement->Type->TypeId;
174 gross 1028 numDXNodesPerElement=0;
175     numCells = elements->numElements;
176 gross 2748 if (TypeId==Line2 || TypeId==Line3 || TypeId==Line4 || TypeId==Line3Macro ) {
177 jgs 153 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 gross 2748 } else if (TypeId==Rec4 || TypeId==Rec8 || TypeId==Rec9 || TypeId==Rec12 || TypeId==Rec16 || TypeId==Rec9Macro ) {
185 jgs 153 numDXNodesPerElement = 4;
186     resortIndex=resort[3];
187     strcpy(elemTypeStr, "quads");
188 gross 2748 } else if (TypeId==Tet4 || TypeId==Tet10 || TypeId==Tet16 || TypeId==Tet10Macro ) {
189 jgs 153 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 gross 2748 sprintf(error_msg,"saveDX: Element type %s is not supported by DX",elements->referenceElementSet->referenceElement->Type->Name);
198 jfenwick 3086 Dudley_setError(VALUE_ERROR,error_msg);
199 gross 1028 MEMFREE(isCellCentered);
200     fclose(fileHandle_p);
201 jgs 153 return;
202     }
203    
204 jgs 82 /* positions */
205 ksteube 1312 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 jgs 82 fprintf(fileHandle_p, "\n");
210     }
211 jgs 153 /* connection table */
212 ksteube 1312 NN=elements->numNodes;
213 jgs 153 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 ksteube 1312 for (j = 0; j < numDXNodesPerElement; j++)
216     fprintf(fileHandle_p," %d",mesh_p->Nodes->reducedNodesMapping->target[elements->Nodes[INDEX2(resortIndex[j], i, NN)]]);
217 jgs 153 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 jgs 82
222     /* data */
223 gross 1028 object_count=2;
224 jgs 153 for (i_data =0 ;i_data<num_data;++i_data) {
225     if (! isEmpty(data_pp[i_data])) {
226     object_count++;
227 gross 1028 rank=getDataPointRank(data_pp[i_data]);
228     nComp=getDataPointSize(data_pp[i_data]);
229 jgs 153 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 jfenwick 3086 if (Dudley_Assemble_reducedIntegrationOrder(data_pp[i_data])) {
236 gross 2748 numPointsPerSample=elements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
237 gross 1062 } else {
238 gross 2748 numPointsPerSample=elements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
239 gross 1062 }
240 jgs 153 if (numPointsPerSample>0) {
241     fprintf(fileHandle_p, "items %d data follows\n", numCells);
242     for (i=0;i<elements->numElements;i++) {
243 jfenwick 2770 values=getSampleDataRO(data_pp[i_data],i);
244 jgs 153 for (k=0;k<nComp;k++) {
245 gross 903 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 jgs 153 }
253     fprintf(fileHandle_p, "\n");
254     }
255     fprintf(fileHandle_p, "attribute \"dep\" string \"connections\"\n");
256     }
257     } else {
258 ksteube 1312 fprintf(fileHandle_p, "items %d data follows\n", numPoints);
259     for (i=0;i<numPoints;i++) {
260 jfenwick 2770 values=getSampleDataRO(data_pp[i_data],i);
261 ksteube 1312 for (k=0;k<nComp;k++) fprintf(fileHandle_p, " %g", values[k]);
262     fprintf(fileHandle_p, "\n");
263 jgs 82 }
264 jgs 153 fprintf(fileHandle_p, "attribute \"dep\" string \"positions\"\n");
265 jgs 82 }
266 jgs 153 }
267 jgs 82 }
268    
269     /* and finish it up */
270 jgs 153 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 jgs 82 fprintf(fileHandle_p, "end\n");
287     /* close the file */
288     fclose(fileHandle_p);
289 gross 1028 MEMFREE(isCellCentered);
290 jgs 82 return;
291     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26