/[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 150 by jgs, Thu Sep 15 03:44:45 2005 UTC trunk/dudley/src/Mesh_saveDX.c revision 3259 by jfenwick, Mon Oct 11 01:48:14 2010 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    /*******************************************************
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  /*   Author: Lutz Gross, gross@access.edu.au */  #include "Mesh.h"
22  /*   Version: $Id$ */  #include "Assemble.h"
23    
24  /**************************************************************/  /**************************************************************/
25    
26  #include "Mesh.h"  void Dudley_Mesh_saveDX(const char *filename_p, Dudley_Mesh * mesh_p, const dim_t num_data, char **names_p,
27                escriptDataC * *data_pp)
28    {
29        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        /* if there is no mesh we just return */
49        if (mesh_p == NULL)
50        return;
51        isCellCentered = MEMALLOC(num_data, bool_t);
52        if (Dudley_checkPtr(isCellCentered))
53        return;
54    
55        fileHandle_p = fopen(filename_p, "w");
56        if (fileHandle_p == NULL)
57        {
58        sprintf(error_msg, "File %s could not be opened for writing.", filename_p);
59        MEMFREE(isCellCentered);
60        fclose(fileHandle_p);
61        Dudley_setError(IO_ERROR, error_msg);
62        return;
63        }
64        /* find the mesh type to be written */
65        elementtype = DUDLEY_UNKNOWN;
66        for (i_data = 0; i_data < num_data; ++i_data)
67        {
68        if (!isEmpty(data_pp[i_data]))
69        {
70            switch (getFunctionSpaceType(data_pp[i_data]))
71            {
72            case DUDLEY_REDUCED_NODES:
73            if (elementtype == DUDLEY_UNKNOWN || elementtype == DUDLEY_ELEMENTS)
74            {
75                elementtype = DUDLEY_ELEMENTS;
76            }
77            else
78            {
79                Dudley_setError(TYPE_ERROR, "saveDX: cannot write given data into single file.");
80                MEMFREE(isCellCentered);
81                fclose(fileHandle_p);
82                return;
83            }
84            isCellCentered[i_data] = FALSE;
85            break;
86            case DUDLEY_ELEMENTS:
87            case DUDLEY_REDUCED_ELEMENTS:
88            if (elementtype == DUDLEY_UNKNOWN || elementtype == DUDLEY_ELEMENTS)
89            {
90                elementtype = DUDLEY_ELEMENTS;
91            }
92            else
93            {
94                Dudley_setError(TYPE_ERROR, "saveDX: cannot write given data into single file.");
95                MEMFREE(isCellCentered);
96                fclose(fileHandle_p);
97                return;
98            }
99            isCellCentered[i_data] = TRUE;
100            break;
101            case DUDLEY_FACE_ELEMENTS:
102            case DUDLEY_REDUCED_FACE_ELEMENTS:
103            if (elementtype == DUDLEY_UNKNOWN || elementtype == DUDLEY_FACE_ELEMENTS)
104            {
105                elementtype = DUDLEY_FACE_ELEMENTS;
106            }
107            else
108            {
109                Dudley_setError(TYPE_ERROR, "saveDX: cannot write given data into single file.");
110                MEMFREE(isCellCentered);
111                fclose(fileHandle_p);
112                return;
113            }
114            isCellCentered[i_data] = TRUE;
115            break;
116            case DUDLEY_POINTS:
117            if (elementtype == DUDLEY_UNKNOWN || elementtype == DUDLEY_POINTS)
118            {
119                elementtype = DUDLEY_POINTS;
120            }
121            else
122            {
123                Dudley_setError(TYPE_ERROR, "saveDX: cannot write given data into single file.");
124                MEMFREE(isCellCentered);
125                fclose(fileHandle_p);
126                return;
127            }
128            isCellCentered[i_data] = TRUE;
129            break;
130            default:
131            sprintf(error_msg, "saveDX: I don't know how to handle function space type %d",
132                getFunctionSpaceType(data_pp[i_data]));
133            Dudley_setError(TYPE_ERROR, error_msg);
134            MEMFREE(isCellCentered);
135            fclose(fileHandle_p);
136            return;
137            }
138        }
139        }
140        /* select number of points and the mesh component */
141        numPoints = mesh_p->Nodes->reducedNodesMapping->numTargets;
142        nDim = mesh_p->Nodes->numDim;
143        if (elementtype == DUDLEY_UNKNOWN)
144        elementtype = DUDLEY_ELEMENTS;
145        switch (elementtype)
146        {
147        case DUDLEY_ELEMENTS:
148        elements = mesh_p->Elements;
149        break;
150        case DUDLEY_FACE_ELEMENTS:
151        elements = mesh_p->FaceElements;
152        break;
153        case DUDLEY_POINTS:
154        elements = mesh_p->Points;
155        break;
156        }
157        if (elements == NULL)
158        {
159        Dudley_setError(SYSTEM_ERROR, "saveDX: undefined element file");
160        MEMFREE(isCellCentered);
161        fclose(fileHandle_p);
162        return;
163        }
164    
165  /**************************************************************/      /* map dudley element type to DX element type */
166        TypeId = elements->etype;
167        numDXNodesPerElement = 0;
168        numCells = elements->numElements;
169        if (TypeId == Dudley_Line2)
170        {
171        numDXNodesPerElement = 2;
172        resortIndex = resort[0];
173        strcpy(elemTypeStr, "lines");
174        }
175        else if (TypeId == Dudley_Tri3)
176        {
177        numDXNodesPerElement = 3;
178        resortIndex = resort[1];
179        strcpy(elemTypeStr, "triangles");
180    /*   } else if (TypeId==Rec4 || TypeId==Rec8 || TypeId==Rec9 || TypeId==Rec12 || TypeId==Rec16 || TypeId==Rec9Macro ) {
181         numDXNodesPerElement = 4;
182         resortIndex=resort[3];
183         strcpy(elemTypeStr, "quads");*/
184        }
185        else if (TypeId == Dudley_Tet4)
186        {
187        numDXNodesPerElement = 4;
188        resortIndex = resort[2];
189        strcpy(elemTypeStr, "tetrahedra");
190    /*   } else if (TypeId==Hex8 || TypeId==Hex20 || TypeId==Hex32 ) {
191         numDXNodesPerElement = 8;
192         resortIndex=resort[4];
193         strcpy(elemTypeStr, "cubes");*/
194        }
195        else
196        {
197        sprintf(error_msg, "saveDX: Element type %s is not supported by DX", elements->ename);
198        Dudley_setError(VALUE_ERROR, error_msg);
199        MEMFREE(isCellCentered);
200        fclose(fileHandle_p);
201        return;
202        }
203    
204  void Finley_Mesh_saveDX(const char * filename_p, Finley_Mesh *mesh_p, escriptDataC* data_p) {      /* positions */
205    char error_msg[LenErrorMsg_MAX];      fprintf(fileHandle_p, "object 1 class array type float rank 1 shape %d items %d data follows\n", nDim, numPoints);
206    /* if there is no mesh we just return */      for (i = 0; i < numPoints; i++)
207    if (mesh_p==NULL) return;      {
208    /* some tables needed for reordering */      p = mesh_p->Nodes->reducedNodesMapping->map[i];
209    int resort[6][9]={      for (j = 0; j < nDim; j++)
210                      {0,1},   /* line */          fprintf(fileHandle_p, " %g", mesh_p->Nodes->Coordinates[INDEX2(j, p, nDim)]);
211                      {0,1,2},  /* triagle */      fprintf(fileHandle_p, "\n");
212                      {0,1,2,3}, /* tetrahedron */      }
213                      {0,3,1,2}, /* quadrilateral */      /* connection table */
214                      {3,0,7,4,2,1,6,5}, /* hexahedron */      NN = elements->numNodes;
215                     };      fprintf(fileHandle_p, "object 2 class array type int rank 1 shape %d items %d data follows\n", numDXNodesPerElement,
216    Finley_ElementFile* elements=NULL;          numCells);
217    char elemTypeStr[32];      for (i = 0; i < numCells; i++)
218    int i,j,k,numDXNodesPerElement,*resortIndex,isCellCentered=TRUE,nodetype=FINLEY_DEGREES_OF_FREEDOM;      {
219    double* values,rtmp;      for (j = 0; j < numDXNodesPerElement; j++)
220    int nDim = mesh_p->Nodes->numDim;          fprintf(fileHandle_p, " %d",
221    /* open the file  and check handel */              mesh_p->Nodes->reducedNodesMapping->target[elements->Nodes[INDEX2(resortIndex[j], i, NN)]]);
222    FILE * fileHandle_p = fopen(filename_p, "w");      fprintf(fileHandle_p, "\n");
223    if (fileHandle_p==NULL) {      }
224           sprintf(error_msg,"File %s could not be opened for writing.",filename_p);      fprintf(fileHandle_p, "attribute \"element type\" string \"%s\"\n", elemTypeStr);
225           Finley_setError(IO_ERROR,error_msg);      fprintf(fileHandle_p, "attribute \"ref\" string \"positions\"\n");
          return;  
   }  
   /* positions */  
   fprintf(fileHandle_p, "object 1 class array type float rank 1 shape %d items %d data follows\n",  
                                                                    nDim, mesh_p->Nodes->reducedNumNodes);  
   for (i = 0; i < mesh_p->Nodes->numNodes; i++) {  
     if (mesh_p->Nodes->toReduced[i]>=0) {  
        fprintf(fileHandle_p, "%g", mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)]);  
        for (j = 1; j < nDim; j++) fprintf(fileHandle_p, " %g",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);  
        fprintf(fileHandle_p, "\n");  
     }  
   }  
   /* connections */  
   /* get a pointer to the relevant mesh component */  
   if (isEmpty(data_p)) {  
       elements=mesh_p->Elements;  
   } else {  
       switch(getFunctionSpaceType(data_p)) {  
        case(FINLEY_DEGREES_OF_FREEDOM):  
           nodetype=FINLEY_DEGREES_OF_FREEDOM;  
           isCellCentered=FALSE;  
           elements=mesh_p->Elements;  
           break;  
        case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):  
           nodetype=FINLEY_REDUCED_DEGREES_OF_FREEDOM;  
           isCellCentered=FALSE;  
           elements=mesh_p->Elements;  
           break;  
        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:  
           sprintf(error_msg,"saveDX:Finley does not know anything about function space type %d",getFunctionSpaceType(data_p));  
           Finley_setError(TYPE_ERROR,error_msg);  
           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 {  
           sprintf(error_msg,"saveDX: Element type %s is not supported by DX",elements->ReferenceElement->Type->Name);  
           Finley_setError(VALUE_ERROR,error_msg);  
           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");  
   
   }  
   /* data */  
   if (!isEmpty(data_p)) {  
       int rank=getDataPointRank(data_p);  
       int nComp=getDataPointSize(data_p);  
       fprintf(fileHandle_p, "object 3 class array type float rank %d ", rank);  
       if (0 < rank) {  
          fprintf(fileHandle_p, "shape ");  
          for (i = 0; i < rank; i++) fprintf(fileHandle_p, "%d ", getDataPointShape(data_p,i));  
       }  
       if (isCellCentered) {  
           int numPointsPerSample=elements->ReferenceElement->numQuadNodes;  
           if (numPointsPerSample>0) {  
              fprintf(fileHandle_p, "items %d data follows\n", elements->numElements);  
              for (i=0;i<elements->numElements;i++) {  
                  values=getSampleData(data_p,i);  
                  for (k=0;k<nComp;k++) {  
                      rtmp=0.;  
                      for (j=0;j<numPointsPerSample;j++) rtmp+=values[INDEX2(k,j,nComp)];  
                      fprintf(fileHandle_p, " %g", rtmp/numPointsPerSample);  
                  }  
              fprintf(fileHandle_p, "\n");  
              }  
              fprintf(fileHandle_p, "attribute \"dep\" string \"connections\"\n");  
          }  
       } 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");  
       }  
   }  
   
   /* and finish it up */  
   fprintf(fileHandle_p, "object 4 class field\n");  
   fprintf(fileHandle_p, "component \"positions\" value 1\n");  
   if (elements!=NULL) fprintf(fileHandle_p, "component \"connections\" value 2\n");  
   if (!isEmpty(data_p)) fprintf(fileHandle_p, "component \"data\" value 3\n");  
   fprintf(fileHandle_p, "end\n");  
   /* close the file */  
   fclose(fileHandle_p);  
   return;  
 }  
226    
227  /*      /* data */
228   * $Log$      object_count = 2;
229   * Revision 1.5  2005/09/15 03:44:23  jgs      for (i_data = 0; i_data < num_data; ++i_data)
230   * Merge of development branch dev-02 back to main trunk on 2005-09-15      {
231   *      if (!isEmpty(data_pp[i_data]))
232   * Revision 1.4.2.3  2005/09/09 08:15:17  gross      {
233   * bugs in vtk and dx writer fixed          object_count++;
234   *          rank = getDataPointRank(data_pp[i_data]);
235   * Revision 1.4.2.2  2005/09/08 08:28:39  gross          nComp = getDataPointSize(data_pp[i_data]);
236   * some cleanup in savevtk          fprintf(fileHandle_p, "object %d class array type float rank %d ", object_count, rank);
237   *          if (0 < rank)
238   * Revision 1.4.2.1  2005/09/07 06:26:20  gross          {
239   * the solver from finley are put into the standalone package paso now          fprintf(fileHandle_p, "shape ");
240   *          for (i = 0; i < rank; i++)
241   * Revision 1.4  2005/07/08 04:07:55  jgs              fprintf(fileHandle_p, "%d ", getDataPointShape(data_pp[i_data], i));
242   * Merge of development branch back to main trunk on 2005-07-08          }
243   *          if (isCellCentered[i_data])
244   * Revision 1.1.1.1.2.5  2005/06/29 02:34:53  gross          {
245   * some changes towards 64 integers in finley          if (Dudley_Assemble_reducedIntegrationOrder(data_pp[i_data]))
246   *          {
247   * Revision 1.1.1.1.2.4  2005/03/03 05:01:12  gross              numPointsPerSample = (elements->numLocalDim == 0) ? 0 : 1;
248   * bug in saveDX fixed          }
249   *          else
250   * Revision 1.1.1.1.2.3  2005/02/17 23:43:06  cochrane          {
251   * Fixed error throwing bug.  Default case of switch statement should have ended              numPointsPerSample = (elements->numLocalDim == 0) ? 0 : (elements->numLocalDim + 1);
252   * with return instead of break, hence errors weren't being thrown (but they now          }
253   * should be).          if (numPointsPerSample > 0)
254   *          {
255   * Revision 1.1.1.1.2.2  2005/02/17 05:53:26  gross              fprintf(fileHandle_p, "items %d data follows\n", numCells);
256   * some bug in saveDX fixed: in fact the bug was in              for (i = 0; i < elements->numElements; i++)
257   * DataC/getDataPointShape              {
258   *              values = getSampleDataRO(data_pp[i_data], i);
259   * Revision 1.1.1.1.2.1  2005/02/17 03:23:01  gross              for (k = 0; k < nComp; k++)
260   * some performance improvements in MVM              {
261   *                  if (isExpanded(data_pp[i_data]))
262   * Revision 1.1.1.1  2004/10/26 06:53:57  jgs                  {
263   * initial import of project esys2                  rtmp = 0.;
264   *                  for (j = 0; j < numPointsPerSample; j++)
265   * Revision 1.1  2004/07/27 08:27:11  gross                      rtmp += values[INDEX2(k, j, nComp)];
266   * Finley: saveDX added: now it is possible to write data on boundary and contact elements                  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.150  
changed lines
  Added in v.3259

  ViewVC Help
Powered by ViewVC 1.1.26