/[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 150 - (hide annotations)
Thu Sep 15 03:44:45 2005 UTC (13 years, 10 months ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Mesh_saveDX.c
File MIME type: text/plain
File size: 10259 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-15

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26