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

  ViewVC Help
Powered by ViewVC 1.1.26