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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3114 - (show annotations)
Fri Aug 27 05:26:25 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: 10602 byte(s)
It doesn't pass all tests but this is major progress

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 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 Dudley_Mesh_saveDX(const char * filename_p, Dudley_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 __const double* values;
41 double rtmp;
42 bool_t *isCellCentered=NULL;
43 Dudley_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 (Dudley_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 Dudley_setError(IO_ERROR,error_msg);
58 return;
59 }
60 /* find the mesh type to be written */
61 elementtype=DUDLEY_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 DUDLEY_REDUCED_NODES:
66 if (elementtype==DUDLEY_UNKNOWN || elementtype==DUDLEY_ELEMENTS) {
67 elementtype=DUDLEY_ELEMENTS;
68 } else {
69 Dudley_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 DUDLEY_ELEMENTS:
77 case DUDLEY_REDUCED_ELEMENTS:
78 if (elementtype==DUDLEY_UNKNOWN || elementtype==DUDLEY_ELEMENTS) {
79 elementtype=DUDLEY_ELEMENTS;
80 } else {
81 Dudley_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 DUDLEY_FACE_ELEMENTS:
89 case DUDLEY_REDUCED_FACE_ELEMENTS:
90 if (elementtype==DUDLEY_UNKNOWN || elementtype==DUDLEY_FACE_ELEMENTS) {
91 elementtype=DUDLEY_FACE_ELEMENTS;
92 } else {
93 Dudley_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 DUDLEY_POINTS:
101 if (elementtype==DUDLEY_UNKNOWN || elementtype==DUDLEY_POINTS) {
102 elementtype=DUDLEY_POINTS;
103 } else {
104 Dudley_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 default:
112 sprintf(error_msg,"saveDX: I don't know how to handle function space type %d",getFunctionSpaceType(data_pp[i_data]));
113 Dudley_setError(TYPE_ERROR,error_msg);
114 MEMFREE(isCellCentered);
115 fclose(fileHandle_p);
116 return;
117 }
118 }
119 }
120 /* select number of points and the mesh component */
121 numPoints = mesh_p->Nodes->reducedNodesMapping->numTargets;
122 nDim = mesh_p->Nodes->numDim;
123 if (elementtype==DUDLEY_UNKNOWN) elementtype=DUDLEY_ELEMENTS;
124 switch(elementtype) {
125 case DUDLEY_ELEMENTS:
126 elements=mesh_p->Elements;
127 break;
128 case DUDLEY_FACE_ELEMENTS:
129 elements=mesh_p->FaceElements;
130 break;
131 case DUDLEY_POINTS:
132 elements=mesh_p->Points;
133 break; }
134 if (elements==NULL) {
135 Dudley_setError(SYSTEM_ERROR,"saveDX: undefined element file");
136 MEMFREE(isCellCentered);
137 fclose(fileHandle_p);
138 return;
139 }
140
141 /* map dudley element type to DX element type */
142 TypeId = elements->referenceElementSet->referenceElement->Type->TypeId;
143 numDXNodesPerElement=0;
144 numCells = elements->numElements;
145 if (TypeId==Line2 || TypeId==Line3 || TypeId==Line4 || TypeId==Line3Macro ) {
146 numDXNodesPerElement=2;
147 resortIndex=resort[0];
148 strcpy(elemTypeStr, "lines");
149 } else if (TypeId==Tri3 || TypeId==Tri6 || TypeId==Tri9 || TypeId==Tri10 ) {
150 numDXNodesPerElement = 3;
151 resortIndex=resort[1];
152 strcpy(elemTypeStr, "triangles");
153 /* } else if (TypeId==Rec4 || TypeId==Rec8 || TypeId==Rec9 || TypeId==Rec12 || TypeId==Rec16 || TypeId==Rec9Macro ) {
154 numDXNodesPerElement = 4;
155 resortIndex=resort[3];
156 strcpy(elemTypeStr, "quads");*/
157 } else if (TypeId==Tet4 || TypeId==Tet10 || TypeId==Tet16 || TypeId==Tet10Macro ) {
158 numDXNodesPerElement = 4;
159 resortIndex=resort[2];
160 strcpy(elemTypeStr, "tetrahedra");
161 /* } else if (TypeId==Hex8 || TypeId==Hex20 || TypeId==Hex32 ) {
162 numDXNodesPerElement = 8;
163 resortIndex=resort[4];
164 strcpy(elemTypeStr, "cubes");*/
165 } else {
166 sprintf(error_msg,"saveDX: Element type %s is not supported by DX",elements->referenceElementSet->referenceElement->Type->Name);
167 Dudley_setError(VALUE_ERROR,error_msg);
168 MEMFREE(isCellCentered);
169 fclose(fileHandle_p);
170 return;
171 }
172
173 /* positions */
174 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 fprintf(fileHandle_p, "\n");
179 }
180 /* connection table */
181 NN=elements->numNodes;
182 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 for (j = 0; j < numDXNodesPerElement; j++)
185 fprintf(fileHandle_p," %d",mesh_p->Nodes->reducedNodesMapping->target[elements->Nodes[INDEX2(resortIndex[j], i, NN)]]);
186 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
191 /* data */
192 object_count=2;
193 for (i_data =0 ;i_data<num_data;++i_data) {
194 if (! isEmpty(data_pp[i_data])) {
195 object_count++;
196 rank=getDataPointRank(data_pp[i_data]);
197 nComp=getDataPointSize(data_pp[i_data]);
198 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 if (Dudley_Assemble_reducedIntegrationOrder(data_pp[i_data])) {
205 numPointsPerSample=elements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
206 } else {
207 numPointsPerSample=elements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
208 }
209 if (numPointsPerSample>0) {
210 fprintf(fileHandle_p, "items %d data follows\n", numCells);
211 for (i=0;i<elements->numElements;i++) {
212 values=getSampleDataRO(data_pp[i_data],i);
213 for (k=0;k<nComp;k++) {
214 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 }
222 fprintf(fileHandle_p, "\n");
223 }
224 fprintf(fileHandle_p, "attribute \"dep\" string \"connections\"\n");
225 }
226 } else {
227 fprintf(fileHandle_p, "items %d data follows\n", numPoints);
228 for (i=0;i<numPoints;i++) {
229 values=getSampleDataRO(data_pp[i_data],i);
230 for (k=0;k<nComp;k++) fprintf(fileHandle_p, " %g", values[k]);
231 fprintf(fileHandle_p, "\n");
232 }
233 fprintf(fileHandle_p, "attribute \"dep\" string \"positions\"\n");
234 }
235 }
236 }
237
238 /* and finish it up */
239 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 fprintf(fileHandle_p, "end\n");
256 /* close the file */
257 fclose(fileHandle_p);
258 MEMFREE(isCellCentered);
259 return;
260 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26