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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26