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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/esys2/finley/src/finleyC/Mesh_saveDX.c revision 115 by jgs, Fri Mar 4 07:12:47 2005 UTC trunk/finley/src/Mesh_saveDX.c revision 903 by gross, Fri Nov 17 01:59:49 2006 UTC
# Line 1  Line 1 
1  /* $Id$ */  /*
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    
# Line 6  Line 17 
17    
18  /**************************************************************/  /**************************************************************/
19    
 /*   Copyrights by ACcESS, Australia 2004 */  
20  /*   Author: Lutz Gross, gross@access.edu.au */  /*   Author: Lutz Gross, gross@access.edu.au */
21    /*   Version: $Id$ */
22    
23  /**************************************************************/  /**************************************************************/
24    
 #include "Finley.h"  
 #include "Common.h"  
25  #include "Mesh.h"  #include "Mesh.h"
 #include "escript/Data/DataC.h"  
26    
27  void Finley_Mesh_saveDX(const char * filename_p, Finley_Mesh *mesh_p, escriptDataC* data_p) {  /**************************************************************/
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];
31    /* if there is no mesh we just return */    /* if there is no mesh we just return */
32    if (mesh_p==NULL) return;    if (mesh_p==NULL) return;
33    
34    /* some tables needed for reordering */    /* some tables needed for reordering */
35    int resort[6][9]={    int resort[6][9]={
36                      {0,1},   /* line */                      {0,1},   /* line */
# Line 27  void Finley_Mesh_saveDX(const char * fil Line 39  void Finley_Mesh_saveDX(const char * fil
39                      {0,3,1,2}, /* quadrilateral */                      {0,3,1,2}, /* quadrilateral */
40                      {3,0,7,4,2,1,6,5}, /* hexahedron */                      {3,0,7,4,2,1,6,5}, /* hexahedron */
41                     };                     };
42    Finley_ElementFile* elements=NULL;    int i,j,k,i_data;
   char elemTypeStr[32];  
   int i,j,k,numDXNodesPerElement,*resortIndex,isCellCentered,nodetype;  
   double* values,rtmp;  
   int nDim = mesh_p->Nodes->numDim;  
43    /* open the file  and check handel */    /* open the file  and check handel */
44    
45    FILE * fileHandle_p = fopen(filename_p, "w");    FILE * fileHandle_p = fopen(filename_p, "w");
46    if (fileHandle_p==NULL) {    if (fileHandle_p==NULL) {
47           Finley_ErrorCode=IO_ERROR;      sprintf(error_msg,"File %s could not be opened for writing.",filename_p);
48           sprintf(Finley_ErrorMsg,"File %s could not be opened for writing.",filename_p);      Finley_setError(IO_ERROR,error_msg);
49           return;      return;
50      }
51      /* find the mesh type to be written */
52      int nodetype=FINLEY_DEGREES_OF_FREEDOM;
53      int elementtype=FINLEY_UNKNOWN;
54      bool_t isCellCentered[num_data];
55      for (i_data=0;i_data<num_data;++i_data) {
56         if (! isEmpty(data_pp[i_data])) {
57            switch(getFunctionSpaceType(data_pp[i_data])) {
58               case FINLEY_DEGREES_OF_FREEDOM:
59                 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
60                 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
61                     elementtype=FINLEY_ELEMENTS;
62                 } else {
63                     Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");
64                     return;
65                 }
66                 isCellCentered[i_data]=FALSE;
67                 break;
68               case FINLEY_REDUCED_DEGREES_OF_FREEDOM:
69                 nodetype = FINLEY_REDUCED_DEGREES_OF_FREEDOM;
70                 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
71                     elementtype=FINLEY_ELEMENTS;
72                 } else {
73                     Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");
74                     return;
75                 }
76                 isCellCentered[i_data]=FALSE;
77                 break;
78               case FINLEY_NODES:
79                 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
80                 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
81                     elementtype=FINLEY_ELEMENTS;
82                 } else {
83                     Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");
84                     return;
85                 }
86                 isCellCentered[i_data]=FALSE;
87                 break;
88               case FINLEY_ELEMENTS:
89                 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
90                 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
91                     elementtype=FINLEY_ELEMENTS;
92                 } else {
93                     Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");
94                     return;
95                 }
96                 isCellCentered[i_data]=TRUE;
97                 break;
98               case FINLEY_FACE_ELEMENTS:
99                 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
100                 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_FACE_ELEMENTS) {
101                     elementtype=FINLEY_FACE_ELEMENTS;
102                 } else {
103                     Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");
104                     return;
105                 }
106                 isCellCentered[i_data]=TRUE;
107                 break;
108               case FINLEY_POINTS:
109                 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
110                 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_POINTS) {
111                     elementtype=FINLEY_POINTS;
112                 } else {
113                     Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");
114                     return;
115                 }
116                 isCellCentered[i_data]=TRUE;
117                 break;
118               case FINLEY_CONTACT_ELEMENTS_1:
119                 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
120                 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {
121                     elementtype=FINLEY_CONTACT_ELEMENTS_1;
122                 } else {
123                     Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");
124                     return;
125                 }
126                 isCellCentered[i_data]=TRUE;
127                 break;
128               case FINLEY_CONTACT_ELEMENTS_2:
129                 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
130                 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {
131                     elementtype=FINLEY_CONTACT_ELEMENTS_1;
132                 } else {
133                     Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");
134                     return;
135                 }
136                 isCellCentered[i_data]=TRUE;
137                 break;
138               default:
139                 sprintf(error_msg,"saveDX: Finley does not know anything about function space type %d",getFunctionSpaceType(data_pp[i_data]));
140                 Finley_setError(TYPE_ERROR,error_msg);
141                 return;
142            }
143         }
144    }    }
145      /* select number of points and the mesh component */
146      int numPoints = mesh_p->Nodes->numNodes;
147      int nDim = mesh_p->Nodes->numDim;
148      if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
149           numPoints = mesh_p->Nodes->reducedNumNodes;
150      } else {
151           numPoints = mesh_p->Nodes->numNodes;
152      }
153      if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;
154      Finley_ElementFile* elements=NULL;
155      switch(elementtype) {
156        case FINLEY_ELEMENTS:
157          elements=mesh_p->Elements;
158          break;
159        case FINLEY_FACE_ELEMENTS:
160          elements=mesh_p->FaceElements;
161          break;
162        case FINLEY_POINTS:
163          elements=mesh_p->Points;
164          break;
165        case FINLEY_CONTACT_ELEMENTS_1:
166          elements=mesh_p->ContactElements;
167          break;
168      }
169      if (elements==NULL) {
170         Finley_setError(SYSTEM_ERROR,"saveDX: undefined element file");
171         return;
172      }
173    
174      /* map finley element type to DX element type */
175      ElementTypeId TypeId = elements->ReferenceElement->Type->TypeId;
176      int *resortIndex=NULL;
177      int numDXNodesPerElement=0;
178      int numCells = elements->numElements;
179      char elemTypeStr[32];
180      if (TypeId==Line2 || TypeId==Line3 || TypeId==Line4 ) {
181         numDXNodesPerElement=2;
182         resortIndex=resort[0];
183         strcpy(elemTypeStr, "lines");
184       } else if (TypeId==Tri3 || TypeId==Tri6 || TypeId==Tri9 || TypeId==Tri10 ) {
185         numDXNodesPerElement = 3;
186         resortIndex=resort[1];
187         strcpy(elemTypeStr, "triangles");
188       } else if (TypeId==Rec4 || TypeId==Rec8 || TypeId==Rec9 || TypeId==Rec12 || TypeId==Rec16 ) {
189         numDXNodesPerElement = 4;
190         resortIndex=resort[3];
191         strcpy(elemTypeStr, "quads");
192       } else if (TypeId==Tet4 || TypeId==Tet10 || TypeId==Tet16 ) {
193         numDXNodesPerElement = 4;
194         resortIndex=resort[2];
195         strcpy(elemTypeStr, "tetrahedra");
196       } else if (TypeId==Hex8 || TypeId==Hex20 || TypeId==Hex32 ) {
197         numDXNodesPerElement = 8;
198         resortIndex=resort[4];
199         strcpy(elemTypeStr, "cubes");
200       } else {
201         sprintf(error_msg,"saveDX: Element type %s is not supported by DX",elements->ReferenceElement->Type->Name);
202         Finley_setError(VALUE_ERROR,error_msg);
203         return;
204       }
205    
206    /* positions */    /* positions */
207    fprintf(fileHandle_p, "object 1 class array type float rank 1 shape %d items %d data follows\n",    fprintf(fileHandle_p, "object 1 class array type float rank 1 shape %d items %d data follows\n",nDim, mesh_p->Nodes->reducedNumNodes);
                                                                    nDim, mesh_p->Nodes->reducedNumNodes);  
208    for (i = 0; i < mesh_p->Nodes->numNodes; i++) {    for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
209      if (mesh_p->Nodes->toReduced[i]>=0) {      if (mesh_p->Nodes->toReduced[i]>=0) {
210         fprintf(fileHandle_p, "%g", mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)]);         for (j = 0; j < nDim; j++) fprintf(fileHandle_p, " %g",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);
        for (j = 1; j < nDim; j++) fprintf(fileHandle_p, " %g",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);  
211         fprintf(fileHandle_p, "\n");         fprintf(fileHandle_p, "\n");
212      }      }
213    }    }
214    /* connections */    /* connection table */
215    /* get a pointer to the relevant mesh component */    int NN=elements->ReferenceElement->Type->numNodes;
216    if (isEmpty(data_p)) {    fprintf(fileHandle_p, "object 2 class array type int rank 1 shape %d items %d data follows\n",numDXNodesPerElement, numCells);
217        elements=mesh_p->Elements;    for (i = 0; i < numCells; i++) {
218    } else {        for (j = 0; j < numDXNodesPerElement; j++) fprintf(fileHandle_p," %d",mesh_p->Nodes->toReduced[elements->Nodes[INDEX2(resortIndex[j], i, NN)]]);
219        switch(getFunctionSpaceType(data_p)) {        fprintf(fileHandle_p, "\n");
220         case(FINLEY_DEGREES_OF_FREEDOM):    }
221            nodetype=FINLEY_DEGREES_OF_FREEDOM;    fprintf(fileHandle_p, "attribute \"element type\" string \"%s\"\n",elemTypeStr);
222            isCellCentered=FALSE;    fprintf(fileHandle_p, "attribute \"ref\" string \"positions\"\n");
           elements=mesh_p->Elements;  
           break;  
        case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):  
           nodetype=FINLEY_REDUCED_DEGREES_OF_FREEDOM;  
           isCellCentered=FALSE;  
           break;  
           elements=mesh_p->Elements;  
        case(FINLEY_NODES):  
           nodetype=FINLEY_NODES;  
           isCellCentered=FALSE;  
           elements=mesh_p->Elements;  
           break;  
        case(FINLEY_ELEMENTS):  
           isCellCentered=TRUE;  
           elements=mesh_p->Elements;  
           break;  
        case(FINLEY_FACE_ELEMENTS):  
           isCellCentered=TRUE;  
           elements=mesh_p->FaceElements;  
           break;  
        case(FINLEY_POINTS):  
           isCellCentered=TRUE;  
           elements=mesh_p->Points;  
           break;  
        case(FINLEY_CONTACT_ELEMENTS_1):  
        case(FINLEY_CONTACT_ELEMENTS_2):  
           isCellCentered=TRUE;  
           elements=mesh_p->ContactElements;  
           break;  
        default:  
           Finley_ErrorCode=TYPE_ERROR;  
           sprintf(Finley_ErrorMsg,"Finley does not know anything about function space type %d",getFunctionSpaceType(data_p));  
           return;  
      }  
   }  
   /* if no element table is present jump over the connection table */  
   if (elements!=NULL) {  
        ElementTypeId TypeId = elements->ReferenceElement->Type->TypeId;  
        if (TypeId==Line2 || TypeId==Line3 || TypeId==Line4 ) {  
           numDXNodesPerElement=2;  
           resortIndex=resort[0];  
           strcpy(elemTypeStr, "lines");  
        } else if (TypeId==Tri3 || TypeId==Tri6 || TypeId==Tri9 || TypeId==Tri10 ) {  
           numDXNodesPerElement = 3;  
           resortIndex=resort[1];  
           strcpy(elemTypeStr, "triangles");  
        } else if (TypeId==Rec4 || TypeId==Rec8 || TypeId==Rec9 || TypeId==Rec12 || TypeId==Rec16 ) {  
           numDXNodesPerElement = 4;  
           resortIndex=resort[3];  
           strcpy(elemTypeStr, "quads");  
         } else if (TypeId==Tet4 || TypeId==Tet10 || TypeId==Tet16 ) {  
           numDXNodesPerElement = 4;  
           resortIndex=resort[2];  
           strcpy(elemTypeStr, "tetrahedra");  
         } else if (TypeId==Hex8 || TypeId==Hex20 || TypeId==Hex32 ) {  
           numDXNodesPerElement = 8;  
           resortIndex=resort[4];  
           strcpy(elemTypeStr, "cubes");  
         } else {  
           Finley_ErrorCode=VALUE_ERROR;  
           sprintf(Finley_ErrorMsg, "Element type %s is not supported by DX",elements->ReferenceElement->Type->Name);  
           return;  
         }  
         int NN=elements->ReferenceElement->Type->numNodes;  
         fprintf(fileHandle_p, "object 2 class array type int rank 1 shape %d items %d data follows\n",numDXNodesPerElement, elements->numElements);  
         for (i = 0; i < elements->numElements; i++) {  
           fprintf(fileHandle_p,"%d",mesh_p->Nodes->toReduced[elements->Nodes[INDEX2(resortIndex[0], i, NN)]]);  
           for (j = 1; j < numDXNodesPerElement; j++) {  
              fprintf(fileHandle_p," %d",mesh_p->Nodes->toReduced[elements->Nodes[INDEX2(resortIndex[j], i, NN)]]);  
           }  
           fprintf(fileHandle_p, "\n");  
         }  
         fprintf(fileHandle_p, "attribute \"element type\" string \"%s\"\n",elemTypeStr);  
         fprintf(fileHandle_p, "attribute \"ref\" string \"positions\"\n");  
223    
   }  
224    /* data */    /* data */
225    if (!isEmpty(data_p)) {    int object_count=2;
226        int rank=getDataPointRank(data_p);    for (i_data =0 ;i_data<num_data;++i_data) {
227        int nComp=getDataPointSize(data_p);        if (! isEmpty(data_pp[i_data])) {
228        fprintf(fileHandle_p, "object 3 class array type float rank %d ", rank);           object_count++;
229        if (0 < rank) {           int rank=getDataPointRank(data_pp[i_data]);
230           fprintf(fileHandle_p, "shape ");           int nComp=getDataPointSize(data_pp[i_data]);
231           for (i = 0; i < rank; i++) fprintf(fileHandle_p, "%d ", getDataPointShape(data_p,i));           double* values,rtmp;
232        }           fprintf(fileHandle_p, "object %d class array type float rank %d ",object_count,rank);
233        if (isCellCentered) {           if (0 < rank) {
234            int numPointsPerSample=elements->ReferenceElement->numQuadNodes;              fprintf(fileHandle_p, "shape ");
235            if (numPointsPerSample>0) {              for (i = 0; i < rank; i++) fprintf(fileHandle_p, "%d ", getDataPointShape(data_pp[i_data],i));
236               fprintf(fileHandle_p, "items %d data follows\n", elements->numElements);           }
237               for (i=0;i<elements->numElements;i++) {           if (isCellCentered[i_data]) {
238                   values=getSampleData(data_p,i);               int numPointsPerSample=elements->ReferenceElement->numQuadNodes;
239                   for (k=0;k<nComp;k++) {               if (numPointsPerSample>0) {
240                       rtmp=0.;                  fprintf(fileHandle_p, "items %d data follows\n", numCells);
241                       for (j=0;j<numPointsPerSample;j++) rtmp+=values[INDEX2(k,j,nComp)];                  for (i=0;i<elements->numElements;i++) {
242                       fprintf(fileHandle_p, " %g", rtmp/numPointsPerSample);                      values=getSampleData(data_pp[i_data],i);
243                        for (k=0;k<nComp;k++) {
244                            if ( isExpanded(data_pp[i_data]) ) {
245                                rtmp=0.;
246                                for (j=0;j<numPointsPerSample;j++) rtmp+=values[INDEX2(k,j,nComp)];
247                                fprintf(fileHandle_p, " %g", rtmp/numPointsPerSample);
248                            } else {
249                                fprintf(fileHandle_p, " %g", values[k]);
250                            }
251                        }
252                    fprintf(fileHandle_p, "\n");
253                    }
254                    fprintf(fileHandle_p, "attribute \"dep\" string \"connections\"\n");
255                }
256             } else {
257                 fprintf(fileHandle_p, "items %d data follows\n", mesh_p->Nodes->reducedNumNodes);
258                 for (i=0;i<mesh_p->Nodes->numNodes;i++) {
259                     if (mesh_p->Nodes->toReduced[i]>=0) {
260                        switch (getFunctionSpaceType(data_pp[i_data])) {
261                           case FINLEY_DEGREES_OF_FREEDOM:
262                              values=getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);
263                              break;
264                           case FINLEY_REDUCED_DEGREES_OF_FREEDOM:
265                              values=getSampleData(data_pp[i_data],mesh_p->Nodes->reducedDegreeOfFreedom[i]);
266                              break;
267                           case FINLEY_NODES:
268                              values=getSampleData(data_pp[i_data],i);
269                              break;
270                        }
271                        for (k=0;k<nComp;k++) fprintf(fileHandle_p, " %g", values[k]);
272                    fprintf(fileHandle_p, "\n");
273                   }                   }
              fprintf(fileHandle_p, "\n");  
274               }               }
275               fprintf(fileHandle_p, "attribute \"dep\" string \"connections\"\n");               fprintf(fileHandle_p, "attribute \"dep\" string \"positions\"\n");
276           }           }
277        } else {       }
           fprintf(fileHandle_p, "items %d data follows\n", mesh_p->Nodes->reducedNumNodes);  
           for (i=0;i<mesh_p->Nodes->numNodes;i++) {  
               if (mesh_p->Nodes->toReduced[i]>=0) {  
                  switch (nodetype) {  
                     case(FINLEY_DEGREES_OF_FREEDOM):  
                        values=getSampleData(data_p,mesh_p->Nodes->degreeOfFreedom[i]);  
                        break;  
                     case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):  
                        values=getSampleData(data_p,mesh_p->Nodes->reducedDegreeOfFreedom[i]);  
                        break;  
                     case(FINLEY_NODES):  
                        values=getSampleData(data_p,i);  
                        break;  
                  }  
                  for (k=0;k<nComp;k++) fprintf(fileHandle_p, " %g", values[k]);  
              fprintf(fileHandle_p, "\n");  
               }  
           }  
           fprintf(fileHandle_p, "attribute \"dep\" string \"positions\"\n");  
       }  
278    }    }
279    
280    /* and finish it up */    /* and finish it up */
281    fprintf(fileHandle_p, "object 4 class field\n");    if (num_data==0) {
282    fprintf(fileHandle_p, "component \"positions\" value 1\n");       fprintf(fileHandle_p, "object %d class field\n",object_count+1);
283    if (elements!=NULL) fprintf(fileHandle_p, "component \"connections\" value 2\n");       fprintf(fileHandle_p, "component \"positions\" value 1\n");
284    if (!isEmpty(data_p)) fprintf(fileHandle_p, "component \"data\" value 3\n");       fprintf(fileHandle_p, "component \"connections\" value 2\n");
285      } else {
286         object_count=2;
287         for (i_data=0; i_data<num_data;++i_data) {
288             if (! isEmpty(data_pp[i_data])) {
289                object_count++;
290                fprintf(fileHandle_p, "object \"%s\" class field\n",names_p[i_data]);
291                fprintf(fileHandle_p, "component \"positions\" value 1\n");
292                fprintf(fileHandle_p, "component \"connections\" value 2\n");
293                fprintf(fileHandle_p, "component \"data\" value %d\n",object_count);
294             }
295         }
296      }
297    fprintf(fileHandle_p, "end\n");    fprintf(fileHandle_p, "end\n");
298    /* close the file */    /* close the file */
299    fclose(fileHandle_p);    fclose(fileHandle_p);

Legend:
Removed from v.115  
changed lines
  Added in v.903

  ViewVC Help
Powered by ViewVC 1.1.26