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

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

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

trunk/finley/src/finley/Mesh_saveDX.c revision 201 by jgs, Wed Nov 23 04:10:21 2005 UTC temp/finley/src/Mesh_saveDX.c revision 1387 by trankine, Fri Jan 11 07:45:26 2008 UTC
# Line 1  Line 1 
 /*  
  ******************************************************************************  
  *                                                                            *  
  *       COPYRIGHT  ACcESS 2003,2004,2005 -  All Rights Reserved              *  
  *                                                                            *  
  * This software is the property of ACcESS. No part of this code              *  
  * may be copied in any form or by any means without the expressed written    *  
  * consent of ACcESS.  Copying, use or modification of this software          *  
  * by any unauthorised person is illegal unless that person has a software    *  
  * license agreement with ACcESS.                                             *  
  *                                                                            *  
  ******************************************************************************  
 */  
1    
2    /* $Id$ */
3    
4  /**************************************************************/  /*******************************************************
5     *
6  /*   writes data and mesh in an opendx file */   *           Copyright 2003-2007 by ACceSS MNRF
7     *       Copyright 2007 by University of Queensland
8     *
9     *                http://esscc.uq.edu.au
10     *        Primary Business: Queensland, Australia
11     *  Licensed under the Open Software License version 3.0
12     *     http://www.opensource.org/licenses/osl-3.0.php
13     *
14     *******************************************************/
15    
16  /**************************************************************/  /**************************************************************/
17    
18  /*   Author: Lutz Gross, gross@access.edu.au */  /*   writes data and mesh in an opendx file                      */
19  /*   Version: $Id$ */  /*   the input data needs to be cell centered or on reducedNodes */
20    
21  /**************************************************************/  /**************************************************************/
22    
23  #include "Mesh.h"  #include "Mesh.h"
24    #include "Assemble.h"
25    
26  /**************************************************************/  /**************************************************************/
27    
28  void Finley_Mesh_saveDX(const char * filename_p, Finley_Mesh *mesh_p, const dim_t num_data,char* *names_p,escriptDataC* *data_pp) {  void Finley_Mesh_saveDX(const char * filename_p, Finley_Mesh *mesh_p, const dim_t num_data,char* *names_p,escriptDataC* *data_pp) {
29    char error_msg[LenErrorMsg_MAX];    char error_msg[LenErrorMsg_MAX], elemTypeStr[32];
   /* if there is no mesh we just return */  
   if (mesh_p==NULL) return;  
   
30    /* some tables needed for reordering */    /* some tables needed for reordering */
31    int resort[6][9]={    int resort[6][9]={
32                      {0,1},   /* line */                      {0,1},   /* line */
# Line 41  void Finley_Mesh_saveDX(const char * fil Line 35  void Finley_Mesh_saveDX(const char * fil
35                      {0,3,1,2}, /* quadrilateral */                      {0,3,1,2}, /* quadrilateral */
36                      {3,0,7,4,2,1,6,5}, /* hexahedron */                      {3,0,7,4,2,1,6,5}, /* hexahedron */
37                     };                     };
38    int i,j,k,i_data;    FILE * fileHandle_p = NULL;
39      int i,j,k,i_data, elementtype, numPoints = 0, nDim, *resortIndex=NULL, p,
40          numDXNodesPerElement=0, numCells, NN, object_count, rank, nComp, numPointsPerSample;
41      double* values,rtmp;
42      bool_t *isCellCentered=NULL;
43      Finley_ElementFile* elements=NULL;
44      ElementTypeId TypeId;
45    /* open the file  and check handel */    /* open the file  and check handel */
46    
47    FILE * fileHandle_p = fopen(filename_p, "w");    /* if there is no mesh we just return */
48      if (mesh_p==NULL) return;
49      isCellCentered=MEMALLOC(num_data, bool_t);
50      if (Finley_checkPtr(isCellCentered)) return;
51    
52      fileHandle_p = fopen(filename_p, "w");
53    if (fileHandle_p==NULL) {    if (fileHandle_p==NULL) {
54      sprintf(error_msg,"File %s could not be opened for writing.",filename_p);      sprintf(error_msg,"File %s could not be opened for writing.",filename_p);
55        MEMFREE(isCellCentered);
56        fclose(fileHandle_p);
57      Finley_setError(IO_ERROR,error_msg);      Finley_setError(IO_ERROR,error_msg);
58      return;      return;
59    }    }
60    /* find the mesh type to be written */    /* find the mesh type to be written */
61    int nodetype=FINLEY_DEGREES_OF_FREEDOM;    elementtype=FINLEY_UNKNOWN;
   int elementtype=FINLEY_UNKNOWN;  
   bool_t isCellCentered[num_data];  
62    for (i_data=0;i_data<num_data;++i_data) {    for (i_data=0;i_data<num_data;++i_data) {
63       if (! isEmpty(data_pp[i_data])) {       if (! isEmpty(data_pp[i_data])) {
64          switch(getFunctionSpaceType(data_pp[i_data])) {          switch(getFunctionSpaceType(data_pp[i_data])) {
65             case FINLEY_DEGREES_OF_FREEDOM:             case FINLEY_REDUCED_NODES:
              nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;  
              if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {  
                  elementtype=FINLEY_ELEMENTS;  
              } else {  
                  Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");  
                  return;  
              }  
              isCellCentered[i_data]=FALSE;  
              break;  
            case FINLEY_REDUCED_DEGREES_OF_FREEDOM:  
              nodetype = FINLEY_REDUCED_DEGREES_OF_FREEDOM;  
66               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
67                   elementtype=FINLEY_ELEMENTS;                   elementtype=FINLEY_ELEMENTS;
68               } else {               } else {
69                   Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");                   Finley_setError(TYPE_ERROR,"saveDX: cannot write given data into single file.");
70                   return;                   MEMFREE(isCellCentered);
71               }                   fclose(fileHandle_p);
              isCellCentered[i_data]=FALSE;  
              break;  
            case FINLEY_NODES:  
              nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;  
              if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {  
                  elementtype=FINLEY_ELEMENTS;  
              } else {  
                  Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");  
72                   return;                   return;
73               }               }
74               isCellCentered[i_data]=FALSE;               isCellCentered[i_data]=FALSE;
75               break;               break;
76             case FINLEY_ELEMENTS:             case FINLEY_ELEMENTS:
77               nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;             case FINLEY_REDUCED_ELEMENTS:
78               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
79                   elementtype=FINLEY_ELEMENTS;                   elementtype=FINLEY_ELEMENTS;
80               } else {               } else {
81                   Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");                   Finley_setError(TYPE_ERROR,"saveDX: cannot write given data into single file.");
82                     MEMFREE(isCellCentered);
83                     fclose(fileHandle_p);
84                   return;                   return;
85               }               }
86               isCellCentered[i_data]=TRUE;               isCellCentered[i_data]=TRUE;
87               break;               break;
88             case FINLEY_FACE_ELEMENTS:             case FINLEY_FACE_ELEMENTS:
89               nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;             case FINLEY_REDUCED_FACE_ELEMENTS:
90               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_FACE_ELEMENTS) {               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_FACE_ELEMENTS) {
91                   elementtype=FINLEY_FACE_ELEMENTS;                   elementtype=FINLEY_FACE_ELEMENTS;
92               } else {               } else {
93                   Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");                   Finley_setError(TYPE_ERROR,"saveDX: cannot write given data into single file.");
94                     MEMFREE(isCellCentered);
95                     fclose(fileHandle_p);
96                   return;                   return;
97               }               }
98               isCellCentered[i_data]=TRUE;               isCellCentered[i_data]=TRUE;
99               break;               break;
100             case FINLEY_POINTS:             case FINLEY_POINTS:
              nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;  
101               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_POINTS) {               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_POINTS) {
102                   elementtype=FINLEY_POINTS;                   elementtype=FINLEY_POINTS;
103               } else {               } else {
104                   Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");                   Finley_setError(TYPE_ERROR,"saveDX: cannot write given data into single file.");
105                     MEMFREE(isCellCentered);
106                     fclose(fileHandle_p);
107                   return;                   return;
108               }               }
109               isCellCentered[i_data]=TRUE;               isCellCentered[i_data]=TRUE;
110               break;               break;
111             case FINLEY_CONTACT_ELEMENTS_1:             case FINLEY_CONTACT_ELEMENTS_1:
112               nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;             case FINLEY_REDUCED_CONTACT_ELEMENTS_1:
113               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {
114                   elementtype=FINLEY_CONTACT_ELEMENTS_1;                   elementtype=FINLEY_CONTACT_ELEMENTS_1;
115               } else {               } else {
116                   Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");                   Finley_setError(TYPE_ERROR,"saveDX: cannot write given data into single file.");
117                     MEMFREE(isCellCentered);
118                     fclose(fileHandle_p);
119                   return;                   return;
120               }               }
121               isCellCentered[i_data]=TRUE;               isCellCentered[i_data]=TRUE;
122               break;               break;
123             case FINLEY_CONTACT_ELEMENTS_2:             case FINLEY_CONTACT_ELEMENTS_2:
124               nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;             case FINLEY_REDUCED_CONTACT_ELEMENTS_2:
125               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {
126                   elementtype=FINLEY_CONTACT_ELEMENTS_1;                   elementtype=FINLEY_CONTACT_ELEMENTS_1;
127               } else {               } else {
128                   Finley_setError(TYPE_ERROR,"saveDX: cannot write given data in single file.");                   Finley_setError(TYPE_ERROR,"saveDX: cannot write given data into single file.");
129                     MEMFREE(isCellCentered);
130                     fclose(fileHandle_p);
131                   return;                   return;
132               }               }
133               isCellCentered[i_data]=TRUE;               isCellCentered[i_data]=TRUE;
134               break;               break;
135             default:             default:
136               sprintf(error_msg,"saveDX: Finley does not know anything about function space type %d",getFunctionSpaceType(data_pp[i_data]));               sprintf(error_msg,"saveDX: I don't know how to handle function space type %d",getFunctionSpaceType(data_pp[i_data]));
137               Finley_setError(TYPE_ERROR,error_msg);               Finley_setError(TYPE_ERROR,error_msg);
138                 MEMFREE(isCellCentered);
139                 fclose(fileHandle_p);
140               return;               return;
141          }          }
142       }       }
143    }    }
144    /* select number of points and the mesh component */    /* select number of points and the mesh component */
145    int numPoints = mesh_p->Nodes->numNodes;    numPoints = mesh_p->Nodes->reducedNodesMapping->numTargets;
146    int nDim = mesh_p->Nodes->numDim;    nDim = mesh_p->Nodes->numDim;
   if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {  
        numPoints = mesh_p->Nodes->reducedNumNodes;  
   } else {  
        numPoints = mesh_p->Nodes->numNodes;  
   }  
147    if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;    if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;
   Finley_ElementFile* elements=NULL;  
148    switch(elementtype) {    switch(elementtype) {
149      case FINLEY_ELEMENTS:      case FINLEY_ELEMENTS:
150        elements=mesh_p->Elements;        elements=mesh_p->Elements;
# Line 164  void Finley_Mesh_saveDX(const char * fil Line 155  void Finley_Mesh_saveDX(const char * fil
155      case FINLEY_POINTS:      case FINLEY_POINTS:
156        elements=mesh_p->Points;        elements=mesh_p->Points;
157        break;        break;
158        case FINLEY_CONTACT_ELEMENTS_2:
159          elements=mesh_p->ContactElements;
160          break;
161      case FINLEY_CONTACT_ELEMENTS_1:      case FINLEY_CONTACT_ELEMENTS_1:
162        elements=mesh_p->ContactElements;        elements=mesh_p->ContactElements;
163        break;        break;
164    }    }
165    if (elements==NULL) {    if (elements==NULL) {
166       Finley_setError(SYSTEM_ERROR,"saveDX: undefined element file");       Finley_setError(SYSTEM_ERROR,"saveDX: undefined element file");
167         MEMFREE(isCellCentered);
168         fclose(fileHandle_p);
169       return;       return;
170    }    }
171    
172    /* map finley element type to DX element type */    /* map finley element type to DX element type */
173    ElementTypeId TypeId = elements->ReferenceElement->Type->TypeId;    TypeId = elements->ReferenceElement->Type->TypeId;
174    int *resortIndex=NULL;    numDXNodesPerElement=0;
175    int numDXNodesPerElement=0;    numCells = elements->numElements;
   int numCells = elements->numElements;  
   char elemTypeStr[32];  
176    if (TypeId==Line2 || TypeId==Line3 || TypeId==Line4 ) {    if (TypeId==Line2 || TypeId==Line3 || TypeId==Line4 ) {
177       numDXNodesPerElement=2;       numDXNodesPerElement=2;
178       resortIndex=resort[0];       resortIndex=resort[0];
# Line 202  void Finley_Mesh_saveDX(const char * fil Line 196  void Finley_Mesh_saveDX(const char * fil
196     } else {     } else {
197       sprintf(error_msg,"saveDX: Element type %s is not supported by DX",elements->ReferenceElement->Type->Name);       sprintf(error_msg,"saveDX: Element type %s is not supported by DX",elements->ReferenceElement->Type->Name);
198       Finley_setError(VALUE_ERROR,error_msg);       Finley_setError(VALUE_ERROR,error_msg);
199         MEMFREE(isCellCentered);
200         fclose(fileHandle_p);
201       return;       return;
202     }     }
203    
204    /* positions */    /* positions */
205    fprintf(fileHandle_p, "object 1 class array type float rank 1 shape %d items %d data follows\n",nDim, mesh_p->Nodes->reducedNumNodes);    fprintf(fileHandle_p, "object 1 class array type float rank 1 shape %d items %d data follows\n",nDim, numPoints);
206    for (i = 0; i < mesh_p->Nodes->numNodes; i++) {    for (i = 0; i < numPoints; i++) {
207      if (mesh_p->Nodes->toReduced[i]>=0) {         p=mesh_p->Nodes->reducedNodesMapping->map[i];
208         for (j = 0; j < nDim; j++) fprintf(fileHandle_p, " %g",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);         for (j = 0; j < nDim; j++) fprintf(fileHandle_p, " %g",mesh_p->Nodes->Coordinates[INDEX2(j, p, nDim)]);
209         fprintf(fileHandle_p, "\n");         fprintf(fileHandle_p, "\n");
     }  
210    }    }
211    /* connection table */    /* connection table */
212    int NN=elements->ReferenceElement->Type->numNodes;    NN=elements->numNodes;
213    fprintf(fileHandle_p, "object 2 class array type int rank 1 shape %d items %d data follows\n",numDXNodesPerElement, numCells);    fprintf(fileHandle_p, "object 2 class array type int rank 1 shape %d items %d data follows\n",numDXNodesPerElement, numCells);
214    for (i = 0; i < numCells; i++) {    for (i = 0; i < numCells; i++) {
215        for (j = 0; j < numDXNodesPerElement; j++) fprintf(fileHandle_p," %d",mesh_p->Nodes->toReduced[elements->Nodes[INDEX2(resortIndex[j], i, NN)]]);        for (j = 0; j < numDXNodesPerElement; j++)
216                fprintf(fileHandle_p," %d",mesh_p->Nodes->reducedNodesMapping->target[elements->Nodes[INDEX2(resortIndex[j], i, NN)]]);
217        fprintf(fileHandle_p, "\n");        fprintf(fileHandle_p, "\n");
218    }    }
219    fprintf(fileHandle_p, "attribute \"element type\" string \"%s\"\n",elemTypeStr);    fprintf(fileHandle_p, "attribute \"element type\" string \"%s\"\n",elemTypeStr);
220    fprintf(fileHandle_p, "attribute \"ref\" string \"positions\"\n");    fprintf(fileHandle_p, "attribute \"ref\" string \"positions\"\n");
221    
222    /* data */    /* data */
223    int object_count=2;    object_count=2;
224    for (i_data =0 ;i_data<num_data;++i_data) {    for (i_data =0 ;i_data<num_data;++i_data) {
225        if (! isEmpty(data_pp[i_data])) {        if (! isEmpty(data_pp[i_data])) {
226           object_count++;           object_count++;
227           int rank=getDataPointRank(data_pp[i_data]);           rank=getDataPointRank(data_pp[i_data]);
228           int nComp=getDataPointSize(data_pp[i_data]);           nComp=getDataPointSize(data_pp[i_data]);
          double* values,rtmp;  
229           fprintf(fileHandle_p, "object %d class array type float rank %d ",object_count,rank);           fprintf(fileHandle_p, "object %d class array type float rank %d ",object_count,rank);
230           if (0 < rank) {           if (0 < rank) {
231              fprintf(fileHandle_p, "shape ");              fprintf(fileHandle_p, "shape ");
232              for (i = 0; i < rank; i++) fprintf(fileHandle_p, "%d ", getDataPointShape(data_pp[i_data],i));              for (i = 0; i < rank; i++) fprintf(fileHandle_p, "%d ", getDataPointShape(data_pp[i_data],i));
233           }           }
234           if (isCellCentered[i_data]) {           if (isCellCentered[i_data]) {
235               int numPointsPerSample=elements->ReferenceElement->numQuadNodes;               if (Finley_Assemble_reducedIntegrationOrder(data_pp[i_data])) {
236                    numPointsPerSample=elements->ReferenceElementReducedOrder->numQuadNodes;
237                 } else {
238                    numPointsPerSample=elements->ReferenceElement->numQuadNodes;
239                 }
240               if (numPointsPerSample>0) {               if (numPointsPerSample>0) {
241                  fprintf(fileHandle_p, "items %d data follows\n", numCells);                  fprintf(fileHandle_p, "items %d data follows\n", numCells);
242                  for (i=0;i<elements->numElements;i++) {                  for (i=0;i<elements->numElements;i++) {
243                      values=getSampleData(data_pp[i_data],i);                      values=getSampleData(data_pp[i_data],i);
244                      for (k=0;k<nComp;k++) {                      for (k=0;k<nComp;k++) {
245                          rtmp=0.;                          if ( isExpanded(data_pp[i_data]) ) {
246                          for (j=0;j<numPointsPerSample;j++) rtmp+=values[INDEX2(k,j,nComp)];                              rtmp=0.;
247                          fprintf(fileHandle_p, " %g", rtmp/numPointsPerSample);                              for (j=0;j<numPointsPerSample;j++) rtmp+=values[INDEX2(k,j,nComp)];
248                                fprintf(fileHandle_p, " %g", rtmp/numPointsPerSample);
249                            } else {
250                                fprintf(fileHandle_p, " %g", values[k]);
251                            }
252                      }                      }
253                  fprintf(fileHandle_p, "\n");                  fprintf(fileHandle_p, "\n");
254                  }                  }
255                  fprintf(fileHandle_p, "attribute \"dep\" string \"connections\"\n");                  fprintf(fileHandle_p, "attribute \"dep\" string \"connections\"\n");
256              }              }
257           } else {           } else {
258               fprintf(fileHandle_p, "items %d data follows\n", mesh_p->Nodes->reducedNumNodes);               fprintf(fileHandle_p, "items %d data follows\n", numPoints);
259               for (i=0;i<mesh_p->Nodes->numNodes;i++) {               for (i=0;i<numPoints;i++) {
260                   if (mesh_p->Nodes->toReduced[i]>=0) {                     values=getSampleData(data_pp[i_data],i);
261                      switch (getFunctionSpaceType(data_pp[i_data])) {                     for (k=0;k<nComp;k++) fprintf(fileHandle_p, " %g", values[k]);
262                         case FINLEY_DEGREES_OF_FREEDOM:                 fprintf(fileHandle_p, "\n");
                           values=getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);  
                           break;  
                        case FINLEY_REDUCED_DEGREES_OF_FREEDOM:  
                           values=getSampleData(data_pp[i_data],mesh_p->Nodes->reducedDegreeOfFreedom[i]);  
                           break;  
                        case FINLEY_NODES:  
                           values=getSampleData(data_pp[i_data],i);  
                           break;  
                     }  
                     for (k=0;k<nComp;k++) fprintf(fileHandle_p, " %g", values[k]);  
                 fprintf(fileHandle_p, "\n");  
                  }  
263               }               }
264               fprintf(fileHandle_p, "attribute \"dep\" string \"positions\"\n");               fprintf(fileHandle_p, "attribute \"dep\" string \"positions\"\n");
265           }           }
# Line 295  void Finley_Mesh_saveDX(const char * fil Line 286  void Finley_Mesh_saveDX(const char * fil
286    fprintf(fileHandle_p, "end\n");    fprintf(fileHandle_p, "end\n");
287    /* close the file */    /* close the file */
288    fclose(fileHandle_p);    fclose(fileHandle_p);
289      MEMFREE(isCellCentered);
290    return;    return;
291  }  }

Legend:
Removed from v.201  
changed lines
  Added in v.1387

  ViewVC Help
Powered by ViewVC 1.1.26