/[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

revision 3221 by jfenwick, Wed Sep 29 01:00:21 2010 UTC revision 3224 by jfenwick, Wed Sep 29 05:19:37 2010 UTC
# 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 Dudley_Mesh_saveDX(const char * filename_p, Dudley_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    Dudley_ElementFile* elements=NULL;      double rtmp;
43    ElementTypeId TypeId;      bool_t *isCellCentered = NULL;
44    /* open the file  and check handel */      Dudley_ElementFile *elements = NULL;
45        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 (Dudley_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      Dudley_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=DUDLEY_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 DUDLEY_REDUCED_NODES:      /* find the mesh type to be written */
65               if (elementtype==DUDLEY_UNKNOWN || elementtype==DUDLEY_ELEMENTS) {      elementtype = DUDLEY_UNKNOWN;
66                   elementtype=DUDLEY_ELEMENTS;      for (i_data = 0; i_data < num_data; ++i_data)
67               } else {      {
68                   Dudley_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 DUDLEY_ELEMENTS:              elementtype = DUDLEY_ELEMENTS;
76             case DUDLEY_REDUCED_ELEMENTS:          }
77               if (elementtype==DUDLEY_UNKNOWN || elementtype==DUDLEY_ELEMENTS) {          else
78                   elementtype=DUDLEY_ELEMENTS;          {
79               } else {              Dudley_setError(TYPE_ERROR, "saveDX: cannot write given data into single file.");
80                   Dudley_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 DUDLEY_FACE_ELEMENTS:          case DUDLEY_REDUCED_ELEMENTS:
88             case DUDLEY_REDUCED_FACE_ELEMENTS:          if (elementtype == DUDLEY_UNKNOWN || elementtype == DUDLEY_ELEMENTS)
89               if (elementtype==DUDLEY_UNKNOWN || elementtype==DUDLEY_FACE_ELEMENTS) {          {
90                   elementtype=DUDLEY_FACE_ELEMENTS;              elementtype = DUDLEY_ELEMENTS;
91               } else {          }
92                   Dudley_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 DUDLEY_POINTS:          isCellCentered[i_data] = TRUE;
100               if (elementtype==DUDLEY_UNKNOWN || elementtype==DUDLEY_POINTS) {          break;
101                   elementtype=DUDLEY_POINTS;          case DUDLEY_FACE_ELEMENTS:
102               } else {          case DUDLEY_REDUCED_FACE_ELEMENTS:
103                   Dudley_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             default:              MEMFREE(isCellCentered);
111               sprintf(error_msg,"saveDX: I don't know how to handle function space type %d",getFunctionSpaceType(data_pp[i_data]));              fclose(fileHandle_p);
112               Dudley_setError(TYPE_ERROR,error_msg);              return;
113               MEMFREE(isCellCentered);          }
114               fclose(fileHandle_p);          isCellCentered[i_data] = TRUE;
115               return;          break;
116          }          case DUDLEY_POINTS:
117       }          if (elementtype == DUDLEY_UNKNOWN || elementtype == DUDLEY_POINTS)
118    }          {
119    /* select number of points and the mesh component */              elementtype = DUDLEY_POINTS;
120    numPoints = mesh_p->Nodes->reducedNodesMapping->numTargets;          }
121    nDim = mesh_p->Nodes->numDim;          else
122    if (elementtype==DUDLEY_UNKNOWN) elementtype=DUDLEY_ELEMENTS;          {
123    switch(elementtype) {              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:      case DUDLEY_ELEMENTS:
148        elements=mesh_p->Elements;      elements = mesh_p->Elements;
149        break;      break;
150      case DUDLEY_FACE_ELEMENTS:      case DUDLEY_FACE_ELEMENTS:
151        elements=mesh_p->FaceElements;      elements = mesh_p->FaceElements;
152        break;      break;
153      case DUDLEY_POINTS:      case DUDLEY_POINTS:
154        elements=mesh_p->Points;      elements = mesh_p->Points;
155        break;  }      break;
156    if (elements==NULL) {      }
157       Dudley_setError(SYSTEM_ERROR,"saveDX: undefined element file");      if (elements == NULL)
158       MEMFREE(isCellCentered);      {
159       fclose(fileHandle_p);      Dudley_setError(SYSTEM_ERROR, "saveDX: undefined element file");
160       return;      MEMFREE(isCellCentered);
161    }      fclose(fileHandle_p);
162        return;
163    /* map dudley element type to DX element type */      }
164    TypeId = elements->etype;  
165    numDXNodesPerElement=0;      /* map dudley element type to DX element type */
166    numCells = elements->numElements;      TypeId = elements->etype;
167    if (TypeId==Line2 ) {      numDXNodesPerElement = 0;
168       numDXNodesPerElement=2;      numCells = elements->numElements;
169       resortIndex=resort[0];      if (TypeId == Line2)
170       strcpy(elemTypeStr, "lines");      {
171     } else if (TypeId==Tri3 ) {      numDXNodesPerElement = 2;
172       numDXNodesPerElement = 3;      resortIndex = resort[0];
173       resortIndex=resort[1];      strcpy(elemTypeStr, "lines");
174       strcpy(elemTypeStr, "triangles");      }
175        else if (TypeId == 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 ) {  /*   } else if (TypeId==Rec4 || TypeId==Rec8 || TypeId==Rec9 || TypeId==Rec12 || TypeId==Rec16 || TypeId==Rec9Macro ) {
181       numDXNodesPerElement = 4;       numDXNodesPerElement = 4;
182       resortIndex=resort[3];       resortIndex=resort[3];
183       strcpy(elemTypeStr, "quads");*/       strcpy(elemTypeStr, "quads");*/
184     } else if (TypeId==Tet4 ) {      }
185       numDXNodesPerElement = 4;      else if (TypeId == Tet4)
186       resortIndex=resort[2];      {
187       strcpy(elemTypeStr, "tetrahedra");      numDXNodesPerElement = 4;
188        resortIndex = resort[2];
189        strcpy(elemTypeStr, "tetrahedra");
190  /*   } else if (TypeId==Hex8 || TypeId==Hex20 || TypeId==Hex32 ) {  /*   } 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->ename);      else
196       Dudley_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 (Dudley_Assemble_reducedIntegrationOrder(data_pp[i_data])) {          object_count++;
234          numPointsPerSample=(elements->numLocalDim==0)?0:1;          rank = getDataPointRank(data_pp[i_data]);
235               } else {          nComp = getDataPointSize(data_pp[i_data]);
236            numPointsPerSample=(elements->numLocalDim==0)?0:(elements->numLocalDim+1);          fprintf(fileHandle_p, "object %d class array type float rank %d ", object_count, rank);
237               }          if (0 < rank)
238               if (numPointsPerSample>0) {          {
239                  fprintf(fileHandle_p, "items %d data follows\n", numCells);          fprintf(fileHandle_p, "shape ");
240                  for (i=0;i<elements->numElements;i++) {          for (i = 0; i < rank; i++)
241                      values=getSampleDataRO(data_pp[i_data],i);              fprintf(fileHandle_p, "%d ", getDataPointShape(data_pp[i_data], i));
242                      for (k=0;k<nComp;k++) {          }
243                          if ( isExpanded(data_pp[i_data]) ) {          if (isCellCentered[i_data])
244                              rtmp=0.;          {
245                              for (j=0;j<numPointsPerSample;j++) rtmp+=values[INDEX2(k,j,nComp)];          if (Dudley_Assemble_reducedIntegrationOrder(data_pp[i_data]))
246                              fprintf(fileHandle_p, " %g", rtmp/numPointsPerSample);          {
247                          } else {              numPointsPerSample = (elements->numLocalDim == 0) ? 0 : 1;
248                              fprintf(fileHandle_p, " %g", values[k]);          }
249                          }          else
250                      }          {
251                  fprintf(fileHandle_p, "\n");              numPointsPerSample = (elements->numLocalDim == 0) ? 0 : (elements->numLocalDim + 1);
252                  }          }
253                  fprintf(fileHandle_p, "attribute \"dep\" string \"connections\"\n");          if (numPointsPerSample > 0)
254              }          {
255           } else {              fprintf(fileHandle_p, "items %d data follows\n", numCells);
256               fprintf(fileHandle_p, "items %d data follows\n", numPoints);              for (i = 0; i < elements->numElements; i++)
257               for (i=0;i<numPoints;i++) {              {
258                     values=getSampleDataRO(data_pp[i_data],i);              values = getSampleDataRO(data_pp[i_data], i);
259                     for (k=0;k<nComp;k++) fprintf(fileHandle_p, " %g", values[k]);              for (k = 0; k < nComp; k++)
260                 fprintf(fileHandle_p, "\n");              {
261               }                  if (isExpanded(data_pp[i_data]))
262               fprintf(fileHandle_p, "attribute \"dep\" string \"positions\"\n");                  {
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    /* and finish it up */                  }
268    if (num_data==0) {                  else
269       fprintf(fileHandle_p, "object %d class field\n",object_count+1);                  {
270       fprintf(fileHandle_p, "component \"positions\" value 1\n");                  fprintf(fileHandle_p, " %g", values[k]);
271       fprintf(fileHandle_p, "component \"connections\" value 2\n");                  }
272    } else {              }
273       object_count=2;              fprintf(fileHandle_p, "\n");
274       for (i_data=0; i_data<num_data;++i_data) {              }
275           if (! isEmpty(data_pp[i_data])) {              fprintf(fileHandle_p, "attribute \"dep\" string \"connections\"\n");
276              object_count++;          }
277              fprintf(fileHandle_p, "object \"%s\" class field\n",names_p[i_data]);          }
278              fprintf(fileHandle_p, "component \"positions\" value 1\n");          else
279              fprintf(fileHandle_p, "component \"connections\" value 2\n");          {
280              fprintf(fileHandle_p, "component \"data\" value %d\n",object_count);          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    fprintf(fileHandle_p, "end\n");              for (k = 0; k < nComp; k++)
285    /* close the file */              fprintf(fileHandle_p, " %g", values[k]);
286    fclose(fileHandle_p);              fprintf(fileHandle_p, "\n");
287    MEMFREE(isCellCentered);          }
288    return;          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.3221  
changed lines
  Added in v.3224

  ViewVC Help
Powered by ViewVC 1.1.26