/[escript]/trunk/esys2/finley/src/finleyC/Mesh_saveDX.c
ViewVC logotype

Contents of /trunk/esys2/finley/src/finleyC/Mesh_saveDX.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 113 - (show annotations)
Mon Feb 28 07:06:33 2005 UTC (14 years, 1 month ago) by jgs
File MIME type: text/plain
File size: 8542 byte(s)
*** empty log message ***

1 /* $Id$ */
2
3 /**************************************************************/
4
5 /* writes data and mesh in an opendx file */
6
7 /**************************************************************/
8
9 /* Copyrights by ACcESS, Australia 2004 */
10 /* Author: Lutz Gross, gross@access.edu.au */
11
12 /**************************************************************/
13
14 #include "Finley.h"
15 #include "Common.h"
16 #include "Mesh.h"
17 #include "escript/Data/DataC.h"
18
19 void Finley_Mesh_saveDX(const char * filename_p, Finley_Mesh *mesh_p, escriptDataC* data_p) {
20 /* if there is no mesh we just return */
21 if (mesh_p==NULL) return;
22 /* some tables needed for reordering */
23 int resort[6][9]={
24 {0,1}, /* line */
25 {0,1,2}, /* triagle */
26 {0,1,2,3}, /* tetrahedron */
27 {0,3,1,2}, /* quadrilateral */
28 {3,0,7,4,2,1,6,5}, /* hexahedron */
29 };
30 Finley_ElementFile* elements=NULL;
31 char elemTypeStr[32];
32 int i,j,k,numDXNodesPerElement,*resortIndex,isCellCentered,nodetype;
33 double* values,rtmp;
34 int nDim = mesh_p->Nodes->numDim;
35 /* open the file and check handel */
36 FILE * fileHandle_p = fopen(filename_p, "w");
37 if (fileHandle_p==NULL) {
38 Finley_ErrorCode=IO_ERROR;
39 sprintf(Finley_ErrorMsg,"File %s could not be opened for writing.",filename_p);
40 return;
41 }
42 /* positions */
43 fprintf(fileHandle_p, "object 1 class array type float rank 1 shape %d items %d data follows\n",
44 nDim, mesh_p->Nodes->reducedNumNodes);
45 for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
46 if (mesh_p->Nodes->toReduced[i]>=0) {
47 fprintf(fileHandle_p, "%g", mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)]);
48 for (j = 1; j < nDim; j++) fprintf(fileHandle_p, " %g",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);
49 fprintf(fileHandle_p, "\n");
50 }
51 }
52 /* connections */
53 /* get a pointer to the relevant mesh component */
54 if (isEmpty(data_p)) {
55 elements=mesh_p->Elements;
56 } else {
57 switch(getFunctionSpaceType(data_p)) {
58 case(FINLEY_DEGREES_OF_FREEDOM):
59 nodetype=FINLEY_DEGREES_OF_FREEDOM;
60 case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):
61 nodetype=FINLEY_REDUCED_DEGREES_OF_FREEDOM;
62 case(FINLEY_NODES):
63 nodetype=FINLEY_NODES;
64 isCellCentered=FALSE;
65 elements=mesh_p->Elements;
66 break;
67 case(FINLEY_ELEMENTS):
68 isCellCentered=TRUE;
69 elements=mesh_p->Elements;
70 break;
71 case(FINLEY_FACE_ELEMENTS):
72 isCellCentered=TRUE;
73 elements=mesh_p->FaceElements;
74 break;
75 case(FINLEY_POINTS):
76 isCellCentered=TRUE;
77 elements=mesh_p->Points;
78 break;
79 case(FINLEY_CONTACT_ELEMENTS_1):
80 case(FINLEY_CONTACT_ELEMENTS_2):
81 isCellCentered=TRUE;
82 elements=mesh_p->ContactElements;
83 break;
84 default:
85 Finley_ErrorCode=TYPE_ERROR;
86 sprintf(Finley_ErrorMsg,"Finley does not know anything about function space type %d",getFunctionSpaceType(data_p));
87 return;
88 }
89 }
90 /* if no element table is present jump over the connection table */
91 if (elements!=NULL) {
92 ElementTypeId TypeId = elements->ReferenceElement->Type->TypeId;
93 if (TypeId==Line2 || TypeId==Line3 || TypeId==Line4 ) {
94 numDXNodesPerElement=2;
95 resortIndex=resort[0];
96 strcpy(elemTypeStr, "lines");
97 } else if (TypeId==Tri3 || TypeId==Tri6 || TypeId==Tri9 || TypeId==Tri10 ) {
98 numDXNodesPerElement = 3;
99 resortIndex=resort[1];
100 strcpy(elemTypeStr, "triangles");
101 } else if (TypeId==Rec4 || TypeId==Rec8 || TypeId==Rec9 || TypeId==Rec12 || TypeId==Rec16 ) {
102 numDXNodesPerElement = 4;
103 resortIndex=resort[3];
104 strcpy(elemTypeStr, "quads");
105 } else if (TypeId==Tet4 || TypeId==Tet10 || TypeId==Tet16 ) {
106 numDXNodesPerElement = 4;
107 resortIndex=resort[2];
108 strcpy(elemTypeStr, "tetrahedra");
109 } else if (TypeId==Hex8 || TypeId==Hex20 || TypeId==Hex32 ) {
110 numDXNodesPerElement = 8;
111 resortIndex=resort[4];
112 strcpy(elemTypeStr, "cubes");
113 } else {
114 Finley_ErrorCode=VALUE_ERROR;
115 sprintf(Finley_ErrorMsg, "Element type %s is not supported by DX",elements->ReferenceElement->Type->Name);
116 return;
117 }
118 int NN=elements->ReferenceElement->Type->numNodes;
119 fprintf(fileHandle_p, "object 2 class array type int rank 1 shape %d items %d data follows\n",numDXNodesPerElement, elements->numElements);
120 for (i = 0; i < elements->numElements; i++) {
121 fprintf(fileHandle_p,"%d",mesh_p->Nodes->toReduced[mesh_p->Elements->Nodes[INDEX2(resortIndex[0], i, NN)]]);
122 for (j = 1; j < numDXNodesPerElement; j++) {
123 fprintf(fileHandle_p," %d",mesh_p->Nodes->toReduced[mesh_p->Elements->Nodes[INDEX2(resortIndex[j], i, NN)]]);
124 }
125 fprintf(fileHandle_p, "\n");
126 }
127 fprintf(fileHandle_p, "attribute \"element type\" string \"%s\"\n",elemTypeStr);
128 fprintf(fileHandle_p, "attribute \"ref\" string \"positions\"\n");
129
130 }
131 /* data */
132 if (!isEmpty(data_p)) {
133 int rank=getDataPointRank(data_p);
134 int nComp=getDataPointSize(data_p);
135 fprintf(fileHandle_p, "object 3 class array type float rank %d ", rank);
136 if (0 < rank) {
137 fprintf(fileHandle_p, "shape ");
138 for (i = 0; i < rank; i++) fprintf(fileHandle_p, "%d ", getDataPointShape(data_p,i));
139 }
140 if (isCellCentered) {
141 int numPointsPerSample=elements->ReferenceElement->numQuadNodes;
142 if (numPointsPerSample>0) {
143 fprintf(fileHandle_p, "items %d data follows\n", elements->numElements);
144 for (i=0;i<elements->numElements;i++) {
145 values=getSampleData(data_p,i);
146 for (k=0;k<nComp;k++) {
147 rtmp=0.;
148 for (j=0;j<numPointsPerSample;j++) rtmp+=values[INDEX2(k,j,nComp)];
149 fprintf(fileHandle_p, " %g", rtmp/numPointsPerSample);
150 }
151 fprintf(fileHandle_p, "\n");
152 }
153 fprintf(fileHandle_p, "attribute \"dep\" string \"connections\"\n");
154 }
155 } else {
156 fprintf(fileHandle_p, "items %d data follows\n", mesh_p->Nodes->reducedNumNodes);
157 for (i=0;i<mesh_p->Nodes->numNodes;i++) {
158 if (mesh_p->Nodes->toReduced[i]>=0) {
159 switch (nodetype) {
160 case(FINLEY_DEGREES_OF_FREEDOM):
161 values=getSampleData(data_p,mesh_p->Nodes->degreeOfFreedom[i]);
162 break;
163 case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):
164 values=getSampleData(data_p,mesh_p->Nodes->reducedDegreeOfFreedom[i]);
165 break;
166 case(FINLEY_NODES):
167 values=getSampleData(data_p,i);
168 break;
169 }
170 }
171 for (k=0;k<nComp;k++) fprintf(fileHandle_p, " %g", values[k]);
172 fprintf(fileHandle_p, "\n");
173 }
174 fprintf(fileHandle_p, "attribute \"dep\" string \"positions\"\n");
175 }
176 }
177
178 /* and finish it up */
179 fprintf(fileHandle_p, "object 4 class field\n");
180 fprintf(fileHandle_p, "component \"positions\" value 1\n");
181 if (elements!=NULL) fprintf(fileHandle_p, "component \"connections\" value 2\n");
182 if (!isEmpty(data_p)) fprintf(fileHandle_p, "component \"data\" value 3\n");
183 fprintf(fileHandle_p, "end\n");
184 /* close the file */
185 fclose(fileHandle_p);
186 return;
187 }
188
189 /*
190 * $Log$
191 * Revision 1.2 2005/02/28 07:06:33 jgs
192 * *** empty log message ***
193 *
194 * Revision 1.1.1.1.2.3 2005/02/17 23:43:06 cochrane
195 * Fixed error throwing bug. Default case of switch statement should have ended
196 * with return instead of break, hence errors weren't being thrown (but they now
197 * should be).
198 *
199 * Revision 1.1.1.1.2.2 2005/02/17 05:53:26 gross
200 * some bug in saveDX fixed: in fact the bug was in
201 * DataC/getDataPointShape
202 *
203 * Revision 1.1.1.1.2.1 2005/02/17 03:23:01 gross
204 * some performance improvements in MVM
205 *
206 * Revision 1.1.1.1 2004/10/26 06:53:57 jgs
207 * initial import of project esys2
208 *
209 * Revision 1.1 2004/07/27 08:27:11 gross
210 * Finley: saveDX added: now it is possible to write data on boundary and contact elements
211 *
212 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26