/[escript]/trunk-mpi-branch/finley/src/Mesh_saveDX.c
ViewVC logotype

Contents of /trunk-mpi-branch/finley/src/Mesh_saveDX.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1223 - (show annotations)
Fri Aug 3 02:40:39 2007 UTC (11 years, 8 months ago) by gross
File MIME type: text/plain
File size: 11853 byte(s)
first attemt towards an improved MPI version.  

1 /*
2 ************************************************************
3 * Copyright 2006 by ACcESS MNRF *
4 * *
5 * http://www.access.edu.au *
6 * Primary Business: Queensland, Australia *
7 * Licensed under the Open Software License version 3.0 *
8 * http://www.opensource.org/licenses/osl-3.0.php *
9 * *
10 ************************************************************
11 */
12
13
14 /**************************************************************/
15
16 /* writes data and mesh in an opendx file */
17 /* the input data needs to be cell centered or on reducedNodes */
18
19 /**************************************************************/
20
21 /* Author: Lutz Gross, gross@access.edu.au */
22 /* Version: $Id$ */
23
24 /**************************************************************/
25
26 #include "Mesh.h"
27 #include "Assemble.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], elemTypeStr[32];
33 /* some tables needed for reordering */
34 int resort[6][9]={
35 {0,1}, /* line */
36 {0,1,2}, /* triagle */
37 {0,1,2,3}, /* tetrahedron */
38 {0,3,1,2}, /* quadrilateral */
39 {3,0,7,4,2,1,6,5}, /* hexahedron */
40 };
41 FILE * fileHandle_p = NULL;
42 int i,j,k,i_data, elementtype, numPoints = 0, nDim, *resortIndex=NULL, p,
43 numDXNodesPerElement=0, numCells, NN, object_count, rank, nComp, numPointsPerSample;
44 double* values,rtmp;
45 bool_t *isCellCentered=NULL;
46 Finley_ElementFile* elements=NULL;
47 ElementTypeId TypeId;
48 /* open the file and check handel */
49
50 /* if there is no mesh we just return */
51 if (mesh_p==NULL) return;
52 isCellCentered=MEMALLOC(num_data, bool_t);
53 if (Finley_checkPtr(isCellCentered)) return;
54
55 fileHandle_p = fopen(filename_p, "w");
56 if (fileHandle_p==NULL) {
57 sprintf(error_msg,"File %s could not be opened for writing.",filename_p);
58 MEMFREE(isCellCentered);
59 fclose(fileHandle_p);
60 Finley_setError(IO_ERROR,error_msg);
61 return;
62 }
63 /* find the mesh type to be written */
64 elementtype=FINLEY_UNKNOWN;
65 for (i_data=0;i_data<num_data;++i_data) {
66 if (! isEmpty(data_pp[i_data])) {
67 switch(getFunctionSpaceType(data_pp[i_data])) {
68 case FINLEY_REDUCED_NODES:
69 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
70 elementtype=FINLEY_ELEMENTS;
71 } else {
72 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");
73 MEMFREE(isCellCentered);
74 fclose(fileHandle_p);
75 return;
76 }
77 isCellCentered[i_data]=FALSE;
78 break;
79 case FINLEY_ELEMENTS:
80 case FINLEY_REDUCED_ELEMENTS:
81 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
82 elementtype=FINLEY_ELEMENTS;
83 } else {
84 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");
85 MEMFREE(isCellCentered);
86 fclose(fileHandle_p);
87 return;
88 }
89 isCellCentered[i_data]=TRUE;
90 break;
91 case FINLEY_FACE_ELEMENTS:
92 case FINLEY_REDUCED_FACE_ELEMENTS:
93 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_FACE_ELEMENTS) {
94 elementtype=FINLEY_FACE_ELEMENTS;
95 } else {
96 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");
97 MEMFREE(isCellCentered);
98 fclose(fileHandle_p);
99 return;
100 }
101 isCellCentered[i_data]=TRUE;
102 break;
103 case FINLEY_POINTS:
104 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_POINTS) {
105 elementtype=FINLEY_POINTS;
106 } else {
107 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");
108 MEMFREE(isCellCentered);
109 fclose(fileHandle_p);
110 return;
111 }
112 isCellCentered[i_data]=TRUE;
113 break;
114 case FINLEY_CONTACT_ELEMENTS_1:
115 case FINLEY_REDUCED_CONTACT_ELEMENTS_1:
116 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {
117 elementtype=FINLEY_CONTACT_ELEMENTS_1;
118 } else {
119 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");
120 MEMFREE(isCellCentered);
121 fclose(fileHandle_p);
122 return;
123 }
124 isCellCentered[i_data]=TRUE;
125 break;
126 case FINLEY_CONTACT_ELEMENTS_2:
127 case FINLEY_REDUCED_CONTACT_ELEMENTS_2:
128 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {
129 elementtype=FINLEY_CONTACT_ELEMENTS_1;
130 } else {
131 Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");
132 MEMFREE(isCellCentered);
133 fclose(fileHandle_p);
134 return;
135 }
136 isCellCentered[i_data]=TRUE;
137 break;
138 default:
139 sprintf(error_msg,"saveDX: I don't know how to handel function space type %d",getFunctionSpaceType(data_pp[i_data]));
140 Finley_setError(TYPE_ERROR,error_msg);
141 MEMFREE(isCellCentered);
142 fclose(fileHandle_p);
143 return;
144 }
145 }
146 }
147 /* select number of points and the mesh component */
148 numPoints = mesh_p->Nodes->reducedNodesMapping->numTargets;
149 nDim = mesh_p->Nodes->numDim;
150 if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;
151 switch(elementtype) {
152 case FINLEY_ELEMENTS:
153 elements=mesh_p->Elements;
154 break;
155 case FINLEY_FACE_ELEMENTS:
156 elements=mesh_p->FaceElements;
157 break;
158 case FINLEY_POINTS:
159 elements=mesh_p->Points;
160 break;
161 case FINLEY_CONTACT_ELEMENTS_2:
162 elements=mesh_p->ContactElements;
163 break;
164 case FINLEY_CONTACT_ELEMENTS_1:
165 elements=mesh_p->ContactElements;
166 break;
167 }
168 if (elements==NULL) {
169 Finley_setError(SYSTEM_ERROR,"saveDX: undefined element file");
170 MEMFREE(isCellCentered);
171 fclose(fileHandle_p);
172 return;
173 }
174
175 /* map finley element type to DX element type */
176 TypeId = elements->ReferenceElement->Type->TypeId;
177 numDXNodesPerElement=0;
178 numCells = elements->numElements;
179 if (TypeId==Line2 || TypeId==Line3 || TypeId==Line4 ) {
180 numDXNodesPerElement=2;
181 resortIndex=resort[0];
182 strcpy(elemTypeStr, "lines");
183 } else if (TypeId==Tri3 || TypeId==Tri6 || TypeId==Tri9 || TypeId==Tri10 ) {
184 numDXNodesPerElement = 3;
185 resortIndex=resort[1];
186 strcpy(elemTypeStr, "triangles");
187 } else if (TypeId==Rec4 || TypeId==Rec8 || TypeId==Rec9 || TypeId==Rec12 || TypeId==Rec16 ) {
188 numDXNodesPerElement = 4;
189 resortIndex=resort[3];
190 strcpy(elemTypeStr, "quads");
191 } else if (TypeId==Tet4 || TypeId==Tet10 || TypeId==Tet16 ) {
192 numDXNodesPerElement = 4;
193 resortIndex=resort[2];
194 strcpy(elemTypeStr, "tetrahedra");
195 } else if (TypeId==Hex8 || TypeId==Hex20 || TypeId==Hex32 ) {
196 numDXNodesPerElement = 8;
197 resortIndex=resort[4];
198 strcpy(elemTypeStr, "cubes");
199 } else {
200 sprintf(error_msg,"saveDX: Element type %s is not supported by DX",elements->ReferenceElement->Type->Name);
201 Finley_setError(VALUE_ERROR,error_msg);
202 MEMFREE(isCellCentered);
203 fclose(fileHandle_p);
204 return;
205 }
206
207 /* positions */
208 fprintf(fileHandle_p, "object 1 class array type float rank 1 shape %d items %d data follows\n",nDim, numPoints);
209 for (i = 0; i < numPoints; i++) {
210 p=mesh_p->Nodes->reducedNodesMapping->map[i];
211 for (j = 0; j < nDim; j++) fprintf(fileHandle_p, " %g",mesh_p->Nodes->Coordinates[INDEX2(j, p, nDim)]);
212 fprintf(fileHandle_p, "\n");
213 }
214 /* connection table */
215 NN=elements->numNodes;
216 fprintf(fileHandle_p, "object 2 class array type int rank 1 shape %d items %d data follows\n",numDXNodesPerElement, numCells);
217 for (i = 0; i < numCells; i++) {
218 for (j = 0; j < numDXNodesPerElement; j++)
219 fprintf(fileHandle_p," %d",mesh_p->Nodes->reducedNodesMapping->target[elements->Nodes[INDEX2(resortIndex[j], i, NN)]]);
220 fprintf(fileHandle_p, "\n");
221 }
222 fprintf(fileHandle_p, "attribute \"element type\" string \"%s\"\n",elemTypeStr);
223 fprintf(fileHandle_p, "attribute \"ref\" string \"positions\"\n");
224
225 /* data */
226 object_count=2;
227 for (i_data =0 ;i_data<num_data;++i_data) {
228 if (! isEmpty(data_pp[i_data])) {
229 object_count++;
230 rank=getDataPointRank(data_pp[i_data]);
231 nComp=getDataPointSize(data_pp[i_data]);
232 fprintf(fileHandle_p, "object %d class array type float rank %d ",object_count,rank);
233 if (0 < rank) {
234 fprintf(fileHandle_p, "shape ");
235 for (i = 0; i < rank; i++) fprintf(fileHandle_p, "%d ", getDataPointShape(data_pp[i_data],i));
236 }
237 if (isCellCentered[i_data]) {
238 if (Finley_Assemble_reducedIntegrationOrder(data_pp[i_data])) {
239 numPointsPerSample=elements->ReferenceElementReducedOrder->numQuadNodes;
240 } else {
241 numPointsPerSample=elements->ReferenceElement->numQuadNodes;
242 }
243 if (numPointsPerSample>0) {
244 fprintf(fileHandle_p, "items %d data follows\n", numCells);
245 for (i=0;i<elements->numElements;i++) {
246 values=getSampleData(data_pp[i_data],i);
247 for (k=0;k<nComp;k++) {
248 if ( isExpanded(data_pp[i_data]) ) {
249 rtmp=0.;
250 for (j=0;j<numPointsPerSample;j++) rtmp+=values[INDEX2(k,j,nComp)];
251 fprintf(fileHandle_p, " %g", rtmp/numPointsPerSample);
252 } else {
253 fprintf(fileHandle_p, " %g", values[k]);
254 }
255 }
256 fprintf(fileHandle_p, "\n");
257 }
258 fprintf(fileHandle_p, "attribute \"dep\" string \"connections\"\n");
259 }
260 } else {
261 fprintf(fileHandle_p, "items %d data follows\n", numPoints);
262 for (i=0;i<numPoints;i++) {
263 values=getSampleData(data_pp[i_data],i);
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