/[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/esys2/finley/src/finleyC/Mesh_saveDX.c revision 82 by jgs, Tue Oct 26 06:53:54 2004 UTC trunk/dudley/src/Mesh_saveDX.c revision 3259 by jfenwick, Mon Oct 11 01:48:14 2010 UTC
# Line 1  Line 1 
1  /* $Id$ */  
2    /*******************************************************
3    *
4    * Copyright (c) 2003-2010 by University of Queensland
5    * Earth Systems Science Computational Center (ESSCC)
6    * http://www.uq.edu.au/esscc
7    *
8    * Primary Business: Queensland, Australia
9    * Licensed under the Open Software License version 3.0
10    * http://www.opensource.org/licenses/osl-3.0.php
11    *
12    *******************************************************/
13    
14  /**************************************************************/  /**************************************************************/
15    
16  /*   writes data and mesh in an opendx file */  /*   writes data and mesh in an opendx file                      */
17    /*   the input data needs to be cell centered or on reducedNodes */
18    
19  /**************************************************************/  /**************************************************************/
20    
21  /*   Copyrights by ACcESS, Australia 2004 */  #include "Mesh.h"
22  /*   Author: Lutz Gross, gross@access.edu.au */  #include "Assemble.h"
23    
24  /**************************************************************/  /**************************************************************/
25    
26  #include "Finley.h"  void Dudley_Mesh_saveDX(const char *filename_p, Dudley_Mesh * mesh_p, const dim_t num_data, char **names_p,
27  #include "Common.h"              escriptDataC * *data_pp)
28  #include "Mesh.h"  {
29  #include "escript/Data/DataC.h"      char error_msg[LenErrorMsg_MAX], elemTypeStr[32];
30        /* some tables needed for reordering */
31        int resort[6][9] = {
32        {0, 1},         /* line */
33        {0, 1, 2},      /* triagle */
34        {0, 1, 2, 3},       /* tetrahedron */
35        {0, 3, 1, 2},       /* quadrilateral */
36        {3, 0, 7, 4, 2, 1, 6, 5},   /* hexahedron */
37        };
38        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        __const double *values;
42        double rtmp;
43        bool_t *isCellCentered = NULL;
44        Dudley_ElementFile *elements = NULL;
45        Dudley_ElementTypeId TypeId;
46        /* open the file  and check handel */
47    
48  void Finley_Mesh_saveDX(const char * filename_p, Finley_Mesh *mesh_p, escriptDataC* data_p) {      /* if there is no mesh we just return */
49    /* if there is no mesh we just return */      if (mesh_p == NULL)
50    if (mesh_p==NULL) return;      return;
51    /* some tables needed for reordering */      isCellCentered = MEMALLOC(num_data, bool_t);
52    int resort[6][9]={      if (Dudley_checkPtr(isCellCentered))
53                      {0,1},   /* line */      return;
54                      {0,1,2},  /* triagle */  
55                      {0,1,2,3}, /* tetrahedron */      fileHandle_p = fopen(filename_p, "w");
56                      {0,3,1,2}, /* quadrilateral */      if (fileHandle_p == NULL)
57                      {3,0,7,4,2,1,6,5}, /* hexahedron */      {
58                     };      sprintf(error_msg, "File %s could not be opened for writing.", filename_p);
59    Finley_ElementFile* elements=NULL;      MEMFREE(isCellCentered);
60    char elemTypeStr[32];      fclose(fileHandle_p);
61    int i,j,k,numDXNodesPerElement,*resortIndex,isCellCentered,nodetype;      Dudley_setError(IO_ERROR, error_msg);
62    double* values,rtmp;      return;
63    int nDim = mesh_p->Nodes->numDim;      }
64    /* open the file  and check handel */      /* find the mesh type to be written */
65    FILE * fileHandle_p = fopen(filename_p, "w");      elementtype = DUDLEY_UNKNOWN;
66    if (fileHandle_p==NULL) {      for (i_data = 0; i_data < num_data; ++i_data)
67           Finley_ErrorCode=IO_ERROR;      {
68           sprintf(Finley_ErrorMsg,"File %s could not be opened for writing.",filename_p);      if (!isEmpty(data_pp[i_data]))
69           return;      {
70    }          switch (getFunctionSpaceType(data_pp[i_data]))
71    /* positions */          {
72    fprintf(fileHandle_p, "object 1 class array type float rank 1 shape %d items %d data follows\n",          case DUDLEY_REDUCED_NODES:
73                                                                     nDim, mesh_p->Nodes->reducedNumNodes);          if (elementtype == DUDLEY_UNKNOWN || elementtype == DUDLEY_ELEMENTS)
74    for (i = 0; i < mesh_p->Nodes->numNodes; i++) {          {
75      if (mesh_p->Nodes->toReduced[i]>=0) {              elementtype = DUDLEY_ELEMENTS;
76         fprintf(fileHandle_p, "%e", mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)]);          }
77         for (j = 1; j < nDim; j++) fprintf(fileHandle_p, " %f",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);          else
78         fprintf(fileHandle_p, "\n");          {
79      }              Dudley_setError(TYPE_ERROR, "saveDX: cannot write given data into single file.");
80    }              MEMFREE(isCellCentered);
81    /* connections */              fclose(fileHandle_p);
82    /* get a pointer to the relevant mesh component */              return;
83    if (isEmpty(data_p)) {          }
84        elements=mesh_p->Elements;          isCellCentered[i_data] = FALSE;
85    } else {          break;
86        switch(getFunctionSpaceType(data_p)) {          case DUDLEY_ELEMENTS:
87         case(FINLEY_DEGREES_OF_FREEDOM):          case DUDLEY_REDUCED_ELEMENTS:
88            nodetype=FINLEY_DEGREES_OF_FREEDOM;          if (elementtype == DUDLEY_UNKNOWN || elementtype == DUDLEY_ELEMENTS)
89         case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):          {
90            nodetype=FINLEY_REDUCED_DEGREES_OF_FREEDOM;              elementtype = DUDLEY_ELEMENTS;
91         case(FINLEY_NODES):          }
92            nodetype=FINLEY_NODES;          else
93            isCellCentered=FALSE;          {
94            elements=mesh_p->Elements;              Dudley_setError(TYPE_ERROR, "saveDX: cannot write given data into single file.");
95            break;              MEMFREE(isCellCentered);
96         case(FINLEY_ELEMENTS):              fclose(fileHandle_p);
97            isCellCentered=TRUE;              return;
98            elements=mesh_p->Elements;          }
99            break;          isCellCentered[i_data] = TRUE;
100         case(FINLEY_FACE_ELEMENTS):          break;
101            isCellCentered=TRUE;          case DUDLEY_FACE_ELEMENTS:
102            elements=mesh_p->FaceElements;          case DUDLEY_REDUCED_FACE_ELEMENTS:
103            break;          if (elementtype == DUDLEY_UNKNOWN || elementtype == DUDLEY_FACE_ELEMENTS)
104         case(FINLEY_POINTS):          {
105            isCellCentered=TRUE;              elementtype = DUDLEY_FACE_ELEMENTS;
106            elements=mesh_p->Points;          }
107            break;          else
108         case(FINLEY_CONTACT_ELEMENTS_1):          {
109         case(FINLEY_CONTACT_ELEMENTS_2):              Dudley_setError(TYPE_ERROR, "saveDX: cannot write given data into single file.");
110            isCellCentered=TRUE;              MEMFREE(isCellCentered);
111            elements=mesh_p->ContactElements;              fclose(fileHandle_p);
112            break;              return;
113         default:          }
114            Finley_ErrorCode=TYPE_ERROR;          isCellCentered[i_data] = TRUE;
115            sprintf(Finley_ErrorMsg,"Finley does not know anything about function space type %d",getFunctionSpaceType(data_p));          break;
116            break;          case DUDLEY_POINTS:
117       }          if (elementtype == DUDLEY_UNKNOWN || elementtype == DUDLEY_POINTS)
118    }          {
119    /* if no element table is present jump over the connection table */              elementtype = DUDLEY_POINTS;
120    if (elements!=NULL) {          }
121         ElementTypeId TypeId = elements->ReferenceElement->Type->TypeId;          else
122         if (TypeId==Line2 || TypeId==Line3 || TypeId==Line4 ) {          {
123            numDXNodesPerElement=2;              Dudley_setError(TYPE_ERROR, "saveDX: cannot write given data into single file.");
124            resortIndex=resort[0];              MEMFREE(isCellCentered);
125            strcpy(elemTypeStr, "lines");              fclose(fileHandle_p);
126         } else if (TypeId==Tri3 || TypeId==Tri6 || TypeId==Tri9 || TypeId==Tri10 ) {              return;
127            numDXNodesPerElement = 3;          }
128            resortIndex=resort[1];          isCellCentered[i_data] = TRUE;
129            strcpy(elemTypeStr, "triangles");          break;
130         } else if (TypeId==Rec4 || TypeId==Rec8 || TypeId==Rec9 || TypeId==Rec12 || TypeId==Rec16 ) {          default:
131            numDXNodesPerElement = 4;          sprintf(error_msg, "saveDX: I don't know how to handle function space type %d",
132            resortIndex=resort[3];              getFunctionSpaceType(data_pp[i_data]));
133            strcpy(elemTypeStr, "quads");          Dudley_setError(TYPE_ERROR, error_msg);
134          } else if (TypeId==Tet4 || TypeId==Tet10 || TypeId==Tet16 ) {          MEMFREE(isCellCentered);
135            numDXNodesPerElement = 4;          fclose(fileHandle_p);
136            resortIndex=resort[2];          return;
137            strcpy(elemTypeStr, "tetrahedra");          }
138          } else if (TypeId==Hex8 || TypeId==Hex20 || TypeId==Hex32 ) {      }
139            numDXNodesPerElement = 8;      }
140            resortIndex=resort[4];      /* select number of points and the mesh component */
141            strcpy(elemTypeStr, "cubes");      numPoints = mesh_p->Nodes->reducedNodesMapping->numTargets;
142          } else {      nDim = mesh_p->Nodes->numDim;
143            Finley_ErrorCode=VALUE_ERROR;      if (elementtype == DUDLEY_UNKNOWN)
144            sprintf(Finley_ErrorMsg, "Element type %s is not supported by DX",elements->ReferenceElement->Type->Name);      elementtype = DUDLEY_ELEMENTS;
145            return;      switch (elementtype)
146          }      {
147          int NN=elements->ReferenceElement->Type->numNodes;      case DUDLEY_ELEMENTS:
148          fprintf(fileHandle_p, "object 2 class array type int rank 1 shape %d items %d data follows\n",numDXNodesPerElement, elements->numElements);      elements = mesh_p->Elements;
149          for (i = 0; i < elements->numElements; i++) {      break;
150            fprintf(fileHandle_p,"%d",mesh_p->Nodes->toReduced[mesh_p->Elements->Nodes[INDEX2(resortIndex[0], i, NN)]]);      case DUDLEY_FACE_ELEMENTS:
151            for (j = 1; j < numDXNodesPerElement; j++) {      elements = mesh_p->FaceElements;
152               fprintf(fileHandle_p," %d",mesh_p->Nodes->toReduced[mesh_p->Elements->Nodes[INDEX2(resortIndex[j], i, NN)]]);      break;
153            }      case DUDLEY_POINTS:
154            fprintf(fileHandle_p, "\n");      elements = mesh_p->Points;
155          }      break;
156          fprintf(fileHandle_p, "attribute \"element type\" string \"%s\"\n",elemTypeStr);      }
157          fprintf(fileHandle_p, "attribute \"ref\" string \"positions\"\n");      if (elements == NULL)
158        {
159    }      Dudley_setError(SYSTEM_ERROR, "saveDX: undefined element file");
160    /* data */      MEMFREE(isCellCentered);
161    if (!isEmpty(data_p)) {      fclose(fileHandle_p);
162        int rank=getDataPointRank(data_p);      return;
163        int* shape=getDataPointShape(data_p);      }
164        int nComp=getDataPointSize(data_p);  
165        fprintf(fileHandle_p, "object 3 class array type float rank %d ", rank);      /* map dudley element type to DX element type */
166        if (0 < rank) {      TypeId = elements->etype;
167           fprintf(fileHandle_p, "shape ");      numDXNodesPerElement = 0;
168           for (i = 0; i < rank; i++) fprintf(fileHandle_p, "%d ", shape[i]);      numCells = elements->numElements;
169        }      if (TypeId == Dudley_Line2)
170        if (isCellCentered) {      {
171            int numPointsPerSample=elements->ReferenceElement->numQuadNodes;      numDXNodesPerElement = 2;
172            if (numPointsPerSample) {      resortIndex = resort[0];
173               fprintf(fileHandle_p, "items %d data follows\n", elements->numElements);      strcpy(elemTypeStr, "lines");
174               for (i=0;i<elements->numElements;i++) {      }
175                   values=getSampleData(data_p,i);      else if (TypeId == Dudley_Tri3)
176                   for (k=0;k<nComp;k++) {      {
177                       rtmp=0.;      numDXNodesPerElement = 3;
178                       for (j=0;j<numPointsPerSample;j++) rtmp+=values[INDEX2(k,j,nComp)];      resortIndex = resort[1];
179                       fprintf(fileHandle_p, " %f", rtmp/numPointsPerSample);      strcpy(elemTypeStr, "triangles");
180                   }  /*   } else if (TypeId==Rec4 || TypeId==Rec8 || TypeId==Rec9 || TypeId==Rec12 || TypeId==Rec16 || TypeId==Rec9Macro ) {
181               fprintf(fileHandle_p, "\n");       numDXNodesPerElement = 4;
182               }       resortIndex=resort[3];
183               fprintf(fileHandle_p, "attribute \"dep\" string \"connections\"\n");       strcpy(elemTypeStr, "quads");*/
184           }      }
185        } else {      else if (TypeId == Dudley_Tet4)
186            fprintf(fileHandle_p, "items %d data follows\n", mesh_p->Nodes->reducedNumNodes);      {
187            for (i=0;i<mesh_p->Nodes->numNodes;i++) {      numDXNodesPerElement = 4;
188                if (mesh_p->Nodes->toReduced[i]>=0) {      resortIndex = resort[2];
189                   switch (nodetype) {      strcpy(elemTypeStr, "tetrahedra");
190                      case(FINLEY_DEGREES_OF_FREEDOM):  /*   } else if (TypeId==Hex8 || TypeId==Hex20 || TypeId==Hex32 ) {
191                         values=getSampleData(data_p,mesh_p->Nodes->degreeOfFreedom[i]);       numDXNodesPerElement = 8;
192                         break;       resortIndex=resort[4];
193                      case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):       strcpy(elemTypeStr, "cubes");*/
194                         values=getSampleData(data_p,mesh_p->Nodes->reducedDegreeOfFreedom[i]);      }
195                         break;      else
196                      case(FINLEY_NODES):      {
197                         values=getSampleData(data_p,i);      sprintf(error_msg, "saveDX: Element type %s is not supported by DX", elements->ename);
198                         break;      Dudley_setError(VALUE_ERROR, error_msg);
199                   }      MEMFREE(isCellCentered);
200                }      fclose(fileHandle_p);
201                for (k=0;k<nComp;k++) fprintf(fileHandle_p, " %f", values[k]);      return;
202            fprintf(fileHandle_p, "\n");      }
203            }  
204            fprintf(fileHandle_p, "attribute \"dep\" string \"positions\"\n");      /* positions */
205        }      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 < numPoints; i++)
207        {
208    /* and finish it up */      p = mesh_p->Nodes->reducedNodesMapping->map[i];
209    fprintf(fileHandle_p, "object 4 class field\n");      for (j = 0; j < nDim; j++)
210    fprintf(fileHandle_p, "component \"positions\" value 1\n");          fprintf(fileHandle_p, " %g", mesh_p->Nodes->Coordinates[INDEX2(j, p, nDim)]);
211    if (elements!=NULL) fprintf(fileHandle_p, "component \"connections\" value 2\n");      fprintf(fileHandle_p, "\n");
212    if (!isEmpty(data_p)) fprintf(fileHandle_p, "component \"data\" value 3\n");      }
213    fprintf(fileHandle_p, "end\n");      /* connection table */
214    /* close the file */      NN = elements->numNodes;
215    fclose(fileHandle_p);      fprintf(fileHandle_p, "object 2 class array type int rank 1 shape %d items %d data follows\n", numDXNodesPerElement,
216    return;          numCells);
217  }      for (i = 0; i < numCells; i++)
218        {
219        for (j = 0; j < numDXNodesPerElement; j++)
220            fprintf(fileHandle_p, " %d",
221                mesh_p->Nodes->reducedNodesMapping->target[elements->Nodes[INDEX2(resortIndex[j], i, NN)]]);
222        fprintf(fileHandle_p, "\n");
223        }
224        fprintf(fileHandle_p, "attribute \"element type\" string \"%s\"\n", elemTypeStr);
225        fprintf(fileHandle_p, "attribute \"ref\" string \"positions\"\n");
226    
227  /*      /* data */
228   * $Log$      object_count = 2;
229   * Revision 1.1  2004/10/26 06:53:57  jgs      for (i_data = 0; i_data < num_data; ++i_data)
230   * Initial revision      {
231   *      if (!isEmpty(data_pp[i_data]))
232   * Revision 1.1  2004/07/27 08:27:11  gross      {
233   * Finley: saveDX added: now it is possible to write data on boundary and contact elements          object_count++;
234   *          rank = getDataPointRank(data_pp[i_data]);
235   */          nComp = getDataPointSize(data_pp[i_data]);
236            fprintf(fileHandle_p, "object %d class array type float rank %d ", object_count, rank);
237            if (0 < rank)
238            {
239            fprintf(fileHandle_p, "shape ");
240            for (i = 0; i < rank; i++)
241                fprintf(fileHandle_p, "%d ", getDataPointShape(data_pp[i_data], i));
242            }
243            if (isCellCentered[i_data])
244            {
245            if (Dudley_Assemble_reducedIntegrationOrder(data_pp[i_data]))
246            {
247                numPointsPerSample = (elements->numLocalDim == 0) ? 0 : 1;
248            }
249            else
250            {
251                numPointsPerSample = (elements->numLocalDim == 0) ? 0 : (elements->numLocalDim + 1);
252            }
253            if (numPointsPerSample > 0)
254            {
255                fprintf(fileHandle_p, "items %d data follows\n", numCells);
256                for (i = 0; i < elements->numElements; i++)
257                {
258                values = getSampleDataRO(data_pp[i_data], i);
259                for (k = 0; k < nComp; k++)
260                {
261                    if (isExpanded(data_pp[i_data]))
262                    {
263                    rtmp = 0.;
264                    for (j = 0; j < numPointsPerSample; j++)
265                        rtmp += values[INDEX2(k, j, nComp)];
266                    fprintf(fileHandle_p, " %g", rtmp / numPointsPerSample);
267                    }
268                    else
269                    {
270                    fprintf(fileHandle_p, " %g", values[k]);
271                    }
272                }
273                fprintf(fileHandle_p, "\n");
274                }
275                fprintf(fileHandle_p, "attribute \"dep\" string \"connections\"\n");
276            }
277            }
278            else
279            {
280            fprintf(fileHandle_p, "items %d data follows\n", numPoints);
281            for (i = 0; i < numPoints; i++)
282            {
283                values = getSampleDataRO(data_pp[i_data], i);
284                for (k = 0; k < nComp; k++)
285                fprintf(fileHandle_p, " %g", values[k]);
286                fprintf(fileHandle_p, "\n");
287            }
288            fprintf(fileHandle_p, "attribute \"dep\" string \"positions\"\n");
289            }
290        }
291        }
292    
293        /* and finish it up */
294        if (num_data == 0)
295        {
296        fprintf(fileHandle_p, "object %d class field\n", object_count + 1);
297        fprintf(fileHandle_p, "component \"positions\" value 1\n");
298        fprintf(fileHandle_p, "component \"connections\" value 2\n");
299        }
300        else
301        {
302        object_count = 2;
303        for (i_data = 0; i_data < num_data; ++i_data)
304        {
305            if (!isEmpty(data_pp[i_data]))
306            {
307            object_count++;
308            fprintf(fileHandle_p, "object \"%s\" class field\n", names_p[i_data]);
309            fprintf(fileHandle_p, "component \"positions\" value 1\n");
310            fprintf(fileHandle_p, "component \"connections\" value 2\n");
311            fprintf(fileHandle_p, "component \"data\" value %d\n", object_count);
312            }
313        }
314        }
315        fprintf(fileHandle_p, "end\n");
316        /* close the file */
317        fclose(fileHandle_p);
318        MEMFREE(isCellCentered);
319        return;
320    }

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

  ViewVC Help
Powered by ViewVC 1.1.26