/[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 82 by jgs, Tue Oct 26 06:53:54 2004 UTC trunk/finley/src/Mesh_saveDX.c revision 616 by elspeth, Wed Mar 22 02:46:56 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, "%e", 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, " %f",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         case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):    fprintf(fileHandle_p, "attribute \"ref\" string \"positions\"\n");
           nodetype=FINLEY_REDUCED_DEGREES_OF_FREEDOM;  
        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));  
           break;  
      }  
   }  
   /* 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[mesh_p->Elements->Nodes[INDEX2(resortIndex[0], i, NN)]]);  
           for (j = 1; j < numDXNodesPerElement; j++) {  
              fprintf(fileHandle_p," %d",mesh_p->Nodes->toReduced[mesh_p->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* shape=getDataPointShape(data_p);        if (! isEmpty(data_pp[i_data])) {
228        int nComp=getDataPointSize(data_p);           object_count++;
229        fprintf(fileHandle_p, "object 3 class array type float rank %d ", rank);           int rank=getDataPointRank(data_pp[i_data]);
230        if (0 < rank) {           int nComp=getDataPointSize(data_pp[i_data]);
231           fprintf(fileHandle_p, "shape ");           double* values,rtmp;
232           for (i = 0; i < rank; i++) fprintf(fileHandle_p, "%d ", shape[i]);           fprintf(fileHandle_p, "object %d class array type float rank %d ",object_count,rank);
233        }           if (0 < rank) {
234        if (isCellCentered) {              fprintf(fileHandle_p, "shape ");
235            int numPointsPerSample=elements->ReferenceElement->numQuadNodes;              for (i = 0; i < rank; i++) fprintf(fileHandle_p, "%d ", getDataPointShape(data_pp[i_data],i));
236            if (numPointsPerSample) {           }
237               fprintf(fileHandle_p, "items %d data follows\n", elements->numElements);           if (isCellCentered[i_data]) {
238               for (i=0;i<elements->numElements;i++) {               int numPointsPerSample=elements->ReferenceElement->numQuadNodes;
239                   values=getSampleData(data_p,i);               if (numPointsPerSample>0) {
240                   for (k=0;k<nComp;k++) {                  fprintf(fileHandle_p, "items %d data follows\n", numCells);
241                       rtmp=0.;                  for (i=0;i<elements->numElements;i++) {
242                       for (j=0;j<numPointsPerSample;j++) rtmp+=values[INDEX2(k,j,nComp)];                      values=getSampleData(data_pp[i_data],i);
243                       fprintf(fileHandle_p, " %f", rtmp/numPointsPerSample);                      for (k=0;k<nComp;k++) {
244                            rtmp=0.;
245                            for (j=0;j<numPointsPerSample;j++) rtmp+=values[INDEX2(k,j,nComp)];
246                            fprintf(fileHandle_p, " %g", rtmp/numPointsPerSample);
247                        }
248                    fprintf(fileHandle_p, "\n");
249                    }
250                    fprintf(fileHandle_p, "attribute \"dep\" string \"connections\"\n");
251                }
252             } else {
253                 fprintf(fileHandle_p, "items %d data follows\n", mesh_p->Nodes->reducedNumNodes);
254                 for (i=0;i<mesh_p->Nodes->numNodes;i++) {
255                     if (mesh_p->Nodes->toReduced[i]>=0) {
256                        switch (getFunctionSpaceType(data_pp[i_data])) {
257                           case FINLEY_DEGREES_OF_FREEDOM:
258                              values=getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);
259                              break;
260                           case FINLEY_REDUCED_DEGREES_OF_FREEDOM:
261                              values=getSampleData(data_pp[i_data],mesh_p->Nodes->reducedDegreeOfFreedom[i]);
262                              break;
263                           case FINLEY_NODES:
264                              values=getSampleData(data_pp[i_data],i);
265                              break;
266                        }
267                        for (k=0;k<nComp;k++) fprintf(fileHandle_p, " %g", values[k]);
268                    fprintf(fileHandle_p, "\n");
269                   }                   }
              fprintf(fileHandle_p, "\n");  
270               }               }
271               fprintf(fileHandle_p, "attribute \"dep\" string \"connections\"\n");               fprintf(fileHandle_p, "attribute \"dep\" string \"positions\"\n");
272           }           }
273        } 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, " %f", values[k]);  
           fprintf(fileHandle_p, "\n");  
           }  
           fprintf(fileHandle_p, "attribute \"dep\" string \"positions\"\n");  
       }  
274    }    }
275    
276    /* and finish it up */    /* and finish it up */
277    fprintf(fileHandle_p, "object 4 class field\n");    if (num_data==0) {
278    fprintf(fileHandle_p, "component \"positions\" value 1\n");       fprintf(fileHandle_p, "object %d class field\n",object_count+1);
279    if (elements!=NULL) fprintf(fileHandle_p, "component \"connections\" value 2\n");       fprintf(fileHandle_p, "component \"positions\" value 1\n");
280    if (!isEmpty(data_p)) fprintf(fileHandle_p, "component \"data\" value 3\n");       fprintf(fileHandle_p, "component \"connections\" value 2\n");
281      } else {
282         object_count=2;
283         for (i_data=0; i_data<num_data;++i_data) {
284             if (! isEmpty(data_pp[i_data])) {
285                object_count++;
286                fprintf(fileHandle_p, "object \"%s\" class field\n",names_p[i_data]);
287                fprintf(fileHandle_p, "component \"positions\" value 1\n");
288                fprintf(fileHandle_p, "component \"connections\" value 2\n");
289                fprintf(fileHandle_p, "component \"data\" value %d\n",object_count);
290             }
291         }
292      }
293    fprintf(fileHandle_p, "end\n");    fprintf(fileHandle_p, "end\n");
294    /* close the file */    /* close the file */
295    fclose(fileHandle_p);    fclose(fileHandle_p);
296    return;    return;
297  }  }
   
 /*  
  * $Log$  
  * Revision 1.1  2004/10/26 06:53:57  jgs  
  * Initial revision  
  *  
  * Revision 1.1  2004/07/27 08:27:11  gross  
  * Finley: saveDX added: now it is possible to write data on boundary and contact elements  
  *  
  */  

Legend:
Removed from v.82  
changed lines
  Added in v.616

  ViewVC Help
Powered by ViewVC 1.1.26