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

Annotation of /trunk/finley/src/Mesh_saveDX.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 123 - (hide annotations)
Fri Jul 8 04:08:13 2005 UTC (14 years ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Mesh_saveDX.c
File MIME type: text/plain
File size: 8949 byte(s)
Merge of development branch back to main trunk on 2005-07-08

1 jgs 82 /* $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 jgs 123 int i,j,k,numDXNodesPerElement,*resortIndex,isCellCentered=TRUE,nodetype=FINLEY_DEGREES_OF_FREEDOM;
33 jgs 82 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 jgs 113 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 jgs 82 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 jgs 115 isCellCentered=FALSE;
61     elements=mesh_p->Elements;
62     break;
63 jgs 82 case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):
64     nodetype=FINLEY_REDUCED_DEGREES_OF_FREEDOM;
65 jgs 115 isCellCentered=FALSE;
66     break;
67     elements=mesh_p->Elements;
68 jgs 82 case(FINLEY_NODES):
69     nodetype=FINLEY_NODES;
70     isCellCentered=FALSE;
71     elements=mesh_p->Elements;
72     break;
73     case(FINLEY_ELEMENTS):
74     isCellCentered=TRUE;
75     elements=mesh_p->Elements;
76     break;
77     case(FINLEY_FACE_ELEMENTS):
78     isCellCentered=TRUE;
79     elements=mesh_p->FaceElements;
80     break;
81     case(FINLEY_POINTS):
82     isCellCentered=TRUE;
83     elements=mesh_p->Points;
84     break;
85     case(FINLEY_CONTACT_ELEMENTS_1):
86     case(FINLEY_CONTACT_ELEMENTS_2):
87     isCellCentered=TRUE;
88     elements=mesh_p->ContactElements;
89     break;
90     default:
91     Finley_ErrorCode=TYPE_ERROR;
92     sprintf(Finley_ErrorMsg,"Finley does not know anything about function space type %d",getFunctionSpaceType(data_p));
93 jgs 113 return;
94 jgs 82 }
95     }
96     /* if no element table is present jump over the connection table */
97     if (elements!=NULL) {
98     ElementTypeId TypeId = elements->ReferenceElement->Type->TypeId;
99     if (TypeId==Line2 || TypeId==Line3 || TypeId==Line4 ) {
100     numDXNodesPerElement=2;
101     resortIndex=resort[0];
102     strcpy(elemTypeStr, "lines");
103     } else if (TypeId==Tri3 || TypeId==Tri6 || TypeId==Tri9 || TypeId==Tri10 ) {
104     numDXNodesPerElement = 3;
105     resortIndex=resort[1];
106     strcpy(elemTypeStr, "triangles");
107     } else if (TypeId==Rec4 || TypeId==Rec8 || TypeId==Rec9 || TypeId==Rec12 || TypeId==Rec16 ) {
108     numDXNodesPerElement = 4;
109     resortIndex=resort[3];
110     strcpy(elemTypeStr, "quads");
111     } else if (TypeId==Tet4 || TypeId==Tet10 || TypeId==Tet16 ) {
112     numDXNodesPerElement = 4;
113     resortIndex=resort[2];
114     strcpy(elemTypeStr, "tetrahedra");
115     } else if (TypeId==Hex8 || TypeId==Hex20 || TypeId==Hex32 ) {
116     numDXNodesPerElement = 8;
117     resortIndex=resort[4];
118     strcpy(elemTypeStr, "cubes");
119     } else {
120     Finley_ErrorCode=VALUE_ERROR;
121     sprintf(Finley_ErrorMsg, "Element type %s is not supported by DX",elements->ReferenceElement->Type->Name);
122     return;
123     }
124     int NN=elements->ReferenceElement->Type->numNodes;
125     fprintf(fileHandle_p, "object 2 class array type int rank 1 shape %d items %d data follows\n",numDXNodesPerElement, elements->numElements);
126     for (i = 0; i < elements->numElements; i++) {
127 jgs 115 fprintf(fileHandle_p,"%d",mesh_p->Nodes->toReduced[elements->Nodes[INDEX2(resortIndex[0], i, NN)]]);
128 jgs 82 for (j = 1; j < numDXNodesPerElement; j++) {
129 jgs 115 fprintf(fileHandle_p," %d",mesh_p->Nodes->toReduced[elements->Nodes[INDEX2(resortIndex[j], i, NN)]]);
130 jgs 82 }
131     fprintf(fileHandle_p, "\n");
132     }
133     fprintf(fileHandle_p, "attribute \"element type\" string \"%s\"\n",elemTypeStr);
134     fprintf(fileHandle_p, "attribute \"ref\" string \"positions\"\n");
135    
136     }
137     /* data */
138     if (!isEmpty(data_p)) {
139     int rank=getDataPointRank(data_p);
140     int nComp=getDataPointSize(data_p);
141     fprintf(fileHandle_p, "object 3 class array type float rank %d ", rank);
142     if (0 < rank) {
143     fprintf(fileHandle_p, "shape ");
144 jgs 113 for (i = 0; i < rank; i++) fprintf(fileHandle_p, "%d ", getDataPointShape(data_p,i));
145 jgs 82 }
146     if (isCellCentered) {
147     int numPointsPerSample=elements->ReferenceElement->numQuadNodes;
148 jgs 113 if (numPointsPerSample>0) {
149 jgs 82 fprintf(fileHandle_p, "items %d data follows\n", elements->numElements);
150     for (i=0;i<elements->numElements;i++) {
151     values=getSampleData(data_p,i);
152     for (k=0;k<nComp;k++) {
153     rtmp=0.;
154     for (j=0;j<numPointsPerSample;j++) rtmp+=values[INDEX2(k,j,nComp)];
155 jgs 113 fprintf(fileHandle_p, " %g", rtmp/numPointsPerSample);
156 jgs 82 }
157     fprintf(fileHandle_p, "\n");
158     }
159     fprintf(fileHandle_p, "attribute \"dep\" string \"connections\"\n");
160     }
161     } else {
162     fprintf(fileHandle_p, "items %d data follows\n", mesh_p->Nodes->reducedNumNodes);
163     for (i=0;i<mesh_p->Nodes->numNodes;i++) {
164     if (mesh_p->Nodes->toReduced[i]>=0) {
165     switch (nodetype) {
166     case(FINLEY_DEGREES_OF_FREEDOM):
167     values=getSampleData(data_p,mesh_p->Nodes->degreeOfFreedom[i]);
168     break;
169     case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):
170     values=getSampleData(data_p,mesh_p->Nodes->reducedDegreeOfFreedom[i]);
171     break;
172     case(FINLEY_NODES):
173     values=getSampleData(data_p,i);
174     break;
175     }
176 jgs 115 for (k=0;k<nComp;k++) fprintf(fileHandle_p, " %g", values[k]);
177     fprintf(fileHandle_p, "\n");
178 jgs 82 }
179     }
180     fprintf(fileHandle_p, "attribute \"dep\" string \"positions\"\n");
181     }
182     }
183    
184     /* and finish it up */
185     fprintf(fileHandle_p, "object 4 class field\n");
186     fprintf(fileHandle_p, "component \"positions\" value 1\n");
187     if (elements!=NULL) fprintf(fileHandle_p, "component \"connections\" value 2\n");
188     if (!isEmpty(data_p)) fprintf(fileHandle_p, "component \"data\" value 3\n");
189     fprintf(fileHandle_p, "end\n");
190     /* close the file */
191     fclose(fileHandle_p);
192     return;
193     }
194 jgs 123
195     /*
196     * $Log$
197     * Revision 1.4 2005/07/08 04:07:55 jgs
198     * Merge of development branch back to main trunk on 2005-07-08
199     *
200     * Revision 1.1.1.1.2.5 2005/06/29 02:34:53 gross
201     * some changes towards 64 integers in finley
202     *
203     * Revision 1.1.1.1.2.4 2005/03/03 05:01:12 gross
204     * bug in saveDX fixed
205     *
206     * Revision 1.1.1.1.2.3 2005/02/17 23:43:06 cochrane
207     * Fixed error throwing bug. Default case of switch statement should have ended
208     * with return instead of break, hence errors weren't being thrown (but they now
209     * should be).
210     *
211     * Revision 1.1.1.1.2.2 2005/02/17 05:53:26 gross
212     * some bug in saveDX fixed: in fact the bug was in
213     * DataC/getDataPointShape
214     *
215     * Revision 1.1.1.1.2.1 2005/02/17 03:23:01 gross
216     * some performance improvements in MVM
217     *
218     * Revision 1.1.1.1 2004/10/26 06:53:57 jgs
219     * initial import of project esys2
220     *
221     * Revision 1.1 2004/07/27 08:27:11 gross
222     * Finley: saveDX added: now it is possible to write data on boundary and contact elements
223     *
224     */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26