/[escript]/trunk/esys2/finley/src/finleyC/Mesh_saveVTK.c
ViewVC logotype

Diff of /trunk/esys2/finley/src/finleyC/Mesh_saveVTK.c

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

revision 112 by jgs, Mon Feb 14 04:14:42 2005 UTC revision 113 by jgs, Mon Feb 28 07:06:33 2005 UTC
# Line 15  Line 15 
15  #include "Common.h"  #include "Common.h"
16  #include "Mesh.h"  #include "Mesh.h"
17  #include "escript/Data/DataC.h"  #include "escript/Data/DataC.h"
18    #include "vtkCellType.h"  /* copied from vtk source directory !!! */
19    
20  void Finley_Mesh_saveVTK(const char * filename_p, Finley_Mesh *mesh_p, escriptDataC* data_p) {  void Finley_Mesh_saveVTK(const char * filename_p, Finley_Mesh *mesh_p, escriptDataC* data_p) {
21    /* if there is no mesh we just return */    /* if there is no mesh we just return */
22    if (mesh_p==NULL) return;    if (mesh_p==NULL) return;
   /* some tables needed for reordering */  
   int resort[6][9]={  
                     {0,1},             /* line */  
                     {0,1,2},           /* triangle */  
                     {0,1,2,3},         /* tetrahedron */  
                     {0,3,1,2},         /* quadrilateral */  
                     {3,0,7,4,2,1,6,5}, /* hexahedron */  
                    };  
23    Finley_ElementFile* elements=NULL;    Finley_ElementFile* elements=NULL;
24    char elemTypeStr[32];    char elemTypeStr[32];
25    int i,j,k,numVTKNodesPerElement,*resortIndex,isCellCentered,nodetype;    int i, j, k, numVTKNodesPerElement, isCellCentered, nodetype;
26    double* values,rtmp;    double* values, rtmp;
27    int nDim = mesh_p->Nodes->numDim;    int nDim = mesh_p->Nodes->numDim;
28    /* open the file  and check handle */  
29      /* get a pointer to the relevant mesh component */
30      if (isEmpty(data_p)) {
31        elements=mesh_p->Elements;
32      } else {
33        switch(getFunctionSpaceType(data_p)) {
34        case(FINLEY_DEGREES_OF_FREEDOM):
35          nodetype = FINLEY_DEGREES_OF_FREEDOM;
36          isCellCentered = FALSE;
37          elements = mesh_p->Elements;
38          break;
39        case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):
40          Finley_ErrorCode=VALUE_ERROR;
41          sprintf(Finley_ErrorMsg,
42              "Reduced degrees of freedom is not yet "
43              "implemented for saving vtk files\n");
44          return;
45        case(FINLEY_NODES):
46          nodetype=FINLEY_NODES;
47          isCellCentered=FALSE;
48          elements=mesh_p->Elements;
49          break;
50        case(FINLEY_ELEMENTS):
51          isCellCentered=TRUE;
52          elements=mesh_p->Elements;
53          break;
54        case(FINLEY_FACE_ELEMENTS):
55          isCellCentered=TRUE;
56          elements=mesh_p->FaceElements;
57          break;
58        case(FINLEY_POINTS):
59          isCellCentered=TRUE;
60          elements=mesh_p->Points;
61          break;
62        case(FINLEY_CONTACT_ELEMENTS_1):
63        case(FINLEY_CONTACT_ELEMENTS_2):
64          isCellCentered=TRUE;
65          elements=mesh_p->ContactElements;
66          break;
67        default:
68          Finley_ErrorCode=TYPE_ERROR;
69          sprintf(Finley_ErrorMsg,
70              "Finley does not know anything about function space type %d",
71              getFunctionSpaceType(data_p));
72          return;
73        }
74      }
75    
76      /* the number of points */
77      int numPoints = mesh_p->Nodes->numNodes;
78    
79      /* the number of cells */
80      if (elements == NULL) {
81        Finley_ErrorCode = VALUE_ERROR;
82        sprintf(Finley_ErrorMsg,
83            "elements object is NULL; cannot proceed");
84        return;
85      }
86      int numCells = elements->numElements;  
87      
88      /* open the file and check handle */
89    FILE * fileHandle_p = fopen(filename_p, "w");    FILE * fileHandle_p = fopen(filename_p, "w");
90    if (fileHandle_p==NULL) {    if (fileHandle_p==NULL) {
91           Finley_ErrorCode=IO_ERROR;      Finley_ErrorCode=IO_ERROR;
92           sprintf(Finley_ErrorMsg,"File %s could not be opened for writing.",filename_p);      sprintf(Finley_ErrorMsg,
93           return;          "File %s could not be opened for writing.", filename_p);
94        return;
95    }    }
96    /* xml header */    /* xml header */
97    fprintf(fileHandle_p, "<?xml version=\"1.0\"?>\n");    fprintf(fileHandle_p, "<?xml version=\"1.0\"?>\n");
# Line 48  void Finley_Mesh_saveVTK(const char * fi Line 102  void Finley_Mesh_saveVTK(const char * fi
102    fprintf(fileHandle_p, "<UnstructuredGrid>\n");    fprintf(fileHandle_p, "<UnstructuredGrid>\n");
103    
104    /* is there only one "piece" to the data?? */    /* is there only one "piece" to the data?? */
105    /* fprintf(fileHandle_p,    fprintf(fileHandle_p, "<Piece "
106        "<Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n",);        "NumberOfPoints=\"%d\" "
107        */        "NumberOfCells=\"%d\">\n",
108          numPoints, numCells);
109    /* now for the point data */  
110    /* fprintf(fileHandle_p,    /* now for the points; equivalent to positions section in saveDX() */
111        "<PointData Scalars=\"%s\" Vectors=\"%s\">\n",);    /* "The points element explicitly defines coordinates for each point
112        */     * individually.  It contains one DataArray element describing an array
113       * with three components per value, each specifying the coordinates of one
114    fprintf(fileHandle_p, "</PointData>\n");     * point" - from Vtk User's Guide
115       */
   /* now for the cell data */  
   /* fprintf(fileHandle_p,  
       "<CellData Scalars=\"%s\" Vectors=\"%s\">\n",);  
       */  
   
   fprintf(fileHandle_p, "</CellData>\n");  
   
   /* now for the points */  
116    fprintf(fileHandle_p, "<Points>\n");    fprintf(fileHandle_p, "<Points>\n");
117    fprintf(fileHandle_p, "<DataArray NumberOfComponents=\"%d\">\n");    /*
118       * the reason for this if statement is explained in the long comment below
119       */
120      if (nDim < 3) {
121        fprintf(fileHandle_p, "<DataArray "
122            "NumberOfComponents=\"3\" "
123            "type=\"Float32\" "
124            "format=\"ascii\">\n");
125      } else {
126        fprintf(fileHandle_p, "<DataArray "
127            "NumberOfComponents=\"%d\" "
128            "type=\"Float32\" "
129            "format=\"ascii\">\n",
130            nDim);
131      }
132      for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
133        fprintf(fileHandle_p,
134            "%e", mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)]);
135        for (j = 1; j < nDim; j++) {
136          fprintf(fileHandle_p,
137              " %f",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);
138          /* vtk/mayavi doesn't like 2D data, it likes 3D data with a degenerate
139           * third dimension to handle 2D data (like a sheet of paper).  So, if
140           * nDim is 2, we have to append zeros to the array to get this third
141           * dimension, and keep the visualisers happy.
142           * Indeed, if nDim is less than 3, must pad all empty dimensions, so
143           * that the total number of dims is 3.
144           */
145          if (nDim < 3) {
146        for (k=0; k<3-nDim; k++) {
147          fprintf(fileHandle_p, " 0");
148        }
149          }
150        }
151        fprintf(fileHandle_p, "\n");
152      }
153    fprintf(fileHandle_p, "</DataArray>\n");    fprintf(fileHandle_p, "</DataArray>\n");
154    fprintf(fileHandle_p, "</Points>\n");    fprintf(fileHandle_p, "</Points>\n");
155    
156      /* connections */
157    /* now for the cells */    /* now for the cells */
158    fprintf(fileHandle_p, "<Cells>\n");    /* "The Cells element defines cells explicitly by specifying point
159    /* fprintf(fileHandle_p,     * connectivity and cell types.  It contains three DataArray elements.  The
160        "<DataArray type=\"%s\" Name=\"%s\" format=\"ascii\">\n",);     * first array specifies the point connectivity.  All cells' point lists
161        */     * are concatenated together.  The second array specifies th eoffset into
162    fprintf(fileHandle_p,     * the connectivity array for the end of each cell.  The third array
163        "</DataArray>\n");     * specifies the type of each cell.
164    fprintf(fileHandle_p, "</Cells>\n");     */
165      /* if no element table is present jump over the connection table */
166      int cellType;
167      if (elements!=NULL) {
168        fprintf(fileHandle_p, "<Cells>\n");
169        ElementTypeId TypeId = elements->ReferenceElement->Type->TypeId;
170        switch(TypeId) {
171        case Point1:
172          cellType = VTK_VERTEX;
173          break;
174        case Line2:
175          cellType = VTK_LINE;
176          break;
177        case Line3:
178          cellType = VTK_QUADRATIC_EDGE;
179          break;
180        case Tri3:
181          cellType = VTK_TRIANGLE;
182          break;
183        case Tri6:
184          cellType = VTK_QUADRATIC_TRIANGLE;
185          break;
186        case Rec4:
187          cellType = VTK_QUAD;
188          break;
189        case Rec8:
190          cellType = VTK_QUADRATIC_QUAD;
191          break;
192        case Tet4:
193          cellType = VTK_TETRA;
194          break;
195        case Tet10:
196          cellType = VTK_QUADRATIC_TETRA;
197          break;
198        case Hex8:
199          cellType = VTK_HEXAHEDRON;
200          break;
201        case Hex20:
202          cellType = VTK_QUADRATIC_HEXAHEDRON;
203          break;
204        case Line2Face:
205          cellType = VTK_VERTEX;
206          break;
207        case Line3Face:
208          cellType = VTK_VERTEX;
209          break;
210        case Tri3Face:
211          cellType = VTK_LINE;
212          break;
213        case Tri6Face:
214          cellType = VTK_QUADRATIC_EDGE;
215          break;
216        case Rec4Face:
217          cellType = VTK_LINE;
218          break;
219        case Rec8Face:
220          cellType = VTK_QUADRATIC_EDGE;
221          break;
222        case Tet4Face:
223          cellType = VTK_TRIANGLE;
224          break;
225        case Tet10Face:
226          cellType = VTK_QUADRATIC_TRIANGLE;
227          break;
228        case Hex8Face:
229          cellType = VTK_QUADRATIC_QUAD;
230          break;
231        case Hex20Face:
232          cellType = VTK_QUADRATIC_QUAD;
233          break;
234        case Point1_Contact:
235          cellType = VTK_VERTEX;
236          break;
237        case Line2_Contact:
238          cellType = VTK_LINE;
239          break;
240        case Line3_Contact:
241          cellType = VTK_QUADRATIC_EDGE;
242          break;
243        case Tri3_Contact:
244          cellType = VTK_TRIANGLE;
245          break;
246        case Tri6_Contact:
247          cellType = VTK_QUADRATIC_TRIANGLE;
248          break;
249        case Rec4_Contact:
250          cellType = VTK_QUAD;
251          break;
252        case Rec8_Contact:
253          cellType = VTK_QUADRATIC_QUAD;
254          break;
255        case Line2Face_Contact:
256          cellType = VTK_VERTEX;
257          break;
258        case Line3Face_Contact:
259          cellType = VTK_VERTEX;
260          break;
261        case Tri3Face_Contact:
262          cellType = VTK_LINE;
263          break;
264        case Tri6Face_Contact:
265          cellType = VTK_QUADRATIC_EDGE;
266          break;
267        case Rec4Face_Contact:
268          cellType = VTK_LINE;
269          break;
270        case Rec8Face_Contact:
271          cellType = VTK_QUADRATIC_EDGE;
272          break;
273        case Tet4Face_Contact:
274          cellType = VTK_TRIANGLE;
275          break;
276        case Tet10Face_Contact:
277          cellType = VTK_QUADRATIC_TRIANGLE;
278          break;
279        case Hex8Face_Contact:
280          cellType = VTK_QUAD;
281          break;
282        case Hex20Face_Contact:
283          cellType = VTK_QUADRATIC_QUAD;
284          break;
285        default:
286          Finley_ErrorCode=VALUE_ERROR;
287          sprintf(Finley_ErrorMsg,
288              "Element type %s is not supported by VTK",
289              elements->ReferenceElement->Type->Name);
290          return;
291        }
292    
293        switch(cellType) {
294        case VTK_VERTEX:
295          numVTKNodesPerElement = 1;
296          strcpy(elemTypeStr, "VTK_VERTEX");
297          break;
298        case VTK_LINE:
299          numVTKNodesPerElement = 2;
300          strcpy(elemTypeStr, "VTK_LINE");
301          break;
302        case VTK_TRIANGLE:
303          numVTKNodesPerElement = 3;
304          strcpy(elemTypeStr, "VTK_TRIANGLE");
305          break;
306        case VTK_QUAD:
307          numVTKNodesPerElement = 4;
308          strcpy(elemTypeStr, "VTK_QUAD");
309          break;
310        case VTK_TETRA:
311          numVTKNodesPerElement = 4;
312          strcpy(elemTypeStr, "VTK_TETRA");
313          break;
314        case VTK_HEXAHEDRON:
315          numVTKNodesPerElement = 8;
316          strcpy(elemTypeStr, "VTK_HEXAHEDRON");
317          break;
318        case VTK_QUADRATIC_EDGE:
319          numVTKNodesPerElement = 3;
320          strcpy(elemTypeStr, "VTK_QUADRATIC_EDGE");
321          break;
322        case VTK_QUADRATIC_TRIANGLE:
323          numVTKNodesPerElement = 6;
324          strcpy(elemTypeStr, "VTK_QUADRATIC_TRIANGLE");
325          break;
326        case VTK_QUADRATIC_QUAD:
327          numVTKNodesPerElement = 8;
328          strcpy(elemTypeStr, "VTK_QUADRATIC_QUAD");
329          break;
330        case VTK_QUADRATIC_TETRA:
331          numVTKNodesPerElement = 10;
332          strcpy(elemTypeStr, "VTK_QUADRATIC_TETRA");
333          break;
334        case VTK_QUADRATIC_HEXAHEDRON:
335          numVTKNodesPerElement = 20;
336          strcpy(elemTypeStr, "VTK_QUADRATIC_HEXAHEDRON");
337          break;
338        default:
339          Finley_ErrorCode = VALUE_ERROR;
340          sprintf(Finley_ErrorMsg,
341              "Cell type %d is not supported by VTK", cellType);
342          return;
343        }
344    
345    /* finish off the piece */      /* write out the DataArray element for the connectivity */
346    fprintf(fileHandle_p, "</Piece>\n");      fprintf(fileHandle_p, "<DataArray "
347            "Name=\"connectivity\" "
348            "type=\"Int32\" "
349            "format=\"ascii\">\n");
350        int NN = elements->ReferenceElement->Type->numNodes;
351        for (i = 0; i < numCells; i++) {
352          fprintf(fileHandle_p, "%d ",
353              mesh_p->Elements->Nodes[INDEX2(0, i, NN)]);
354          for (j = 1; j < numVTKNodesPerElement; j++) {
355        fprintf(fileHandle_p,"%d ", mesh_p->Elements->Nodes[INDEX2(j, i, NN)]);
356          }
357        }
358        fprintf(fileHandle_p, "\n");
359        fprintf(fileHandle_p, "</DataArray>\n");
360    
361        /* write out the DataArray element for the offsets */
362        fprintf(fileHandle_p, "<DataArray "
363            "Name=\"offsets\" "
364            "type=\"Int32\" "
365            "format=\"ascii\">\n");
366        for (i=numVTKNodesPerElement; i<=numCells*numVTKNodesPerElement; i+=numVTKNodesPerElement) {
367          fprintf(fileHandle_p, "%d ", i);
368        }
369        fprintf(fileHandle_p, "\n");
370        fprintf(fileHandle_p, "</DataArray>\n");
371    
372    /* positions */      /* write out the DataArray element for the types */
373    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, "<DataArray "
374    for (i = 0; i < mesh_p->Nodes->numNodes; i++) {          "Name=\"types\" "
375      if (mesh_p->Nodes->toReduced[i]>=0) {          "type=\"UInt8\" "
376         fprintf(fileHandle_p, "%e", mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)]);          "format=\"ascii\">\n");
377         for (j = 1; j < nDim; j++) fprintf(fileHandle_p, " %f",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);      for (i=0; i<numCells; i++) {
378         fprintf(fileHandle_p, "\n");        fprintf(fileHandle_p, "%d ", cellType);
379      }      }
380    }      fprintf(fileHandle_p, "\n");
381    /* connections */      fprintf(fileHandle_p, "</DataArray>\n");
   /* 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;  
        case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):  
           nodetype=FINLEY_REDUCED_DEGREES_OF_FREEDOM;  
        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:  
           Finley_ErrorCode=TYPE_ERROR;  
           sprintf(Finley_ErrorMsg,"Finley does not know anything about function space type %d",getFunctionSpaceType(data_p));  
           break;  
      }  
   }  
   /* 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 ) {  
           numVTKNodesPerElement=2;  
           resortIndex=resort[0];  
           strcpy(elemTypeStr, "lines");  
        } else if (TypeId==Tri3 || TypeId==Tri6 || TypeId==Tri9 || TypeId==Tri10 ) {  
           numVTKNodesPerElement = 3;  
           resortIndex=resort[1];  
           strcpy(elemTypeStr, "triangles");  
        } else if (TypeId==Rec4 || TypeId==Rec8 || TypeId==Rec9 || TypeId==Rec12 || TypeId==Rec16 ) {  
           numVTKNodesPerElement = 4;  
           resortIndex=resort[3];  
           strcpy(elemTypeStr, "quads");  
         } else if (TypeId==Tet4 || TypeId==Tet10 || TypeId==Tet16 ) {  
           numVTKNodesPerElement = 4;  
           resortIndex=resort[2];  
           strcpy(elemTypeStr, "tetrahedra");  
         } else if (TypeId==Hex8 || TypeId==Hex20 || TypeId==Hex32 ) {  
           numVTKNodesPerElement = 8;  
           resortIndex=resort[4];  
           strcpy(elemTypeStr, "cubes");  
         } else {  
           Finley_ErrorCode=VALUE_ERROR;  
           sprintf(Finley_ErrorMsg, "Element type %s is not supported by VTK",elements->ReferenceElement->Type->Name);  
           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",numVTKNodesPerElement, elements->numElements);  
         for (i = 0; i < elements->numElements; i++) {  
           fprintf(fileHandle_p,"%d",mesh_p->Nodes->toReduced[mesh_p->Elements->Nodes[INDEX2(resortIndex[0], i, NN)]]);  
           for (j = 1; j < numVTKNodesPerElement; j++) {  
              fprintf(fileHandle_p," %d",mesh_p->Nodes->toReduced[mesh_p->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");  
382    
383        /* finish off the <Cells> element */
384        fprintf(fileHandle_p, "</Cells>\n");
385    }    }
386    
387    /* data */    /* data */
388    if (!isEmpty(data_p)) {    if (!isEmpty(data_p)) {
389        int rank=getDataPointRank(data_p);      int rank = getDataPointRank(data_p);
390        int* shape=getDataPointShape(data_p);      int nComp = getDataPointSize(data_p);
391        int nComp=getDataPointSize(data_p);      /* barf if rank is greater than two */
392        fprintf(fileHandle_p, "object 3 class array type float rank %d ", rank);      if (rank > 2) {
393        if (0 < rank) {        Finley_ErrorCode = VALUE_ERROR;
394           fprintf(fileHandle_p, "shape ");        sprintf(Finley_ErrorMsg,
395           for (i = 0; i < rank; i++) fprintf(fileHandle_p, "%d ", shape[i]);            "Vtk can't handle objects with rank greater than 2\n"
396              "object rank = %d\n", rank);
397          return;
398        }
399        /* if the rank == 0:   --> scalar data
400         * if the rank == 1:   --> vector data
401         * if the rank == 2:   --> tensor data
402         */
403        char dataNameStr[31], dataTypeStr[63];
404        int nCompReqd;   /* the number of components required by vtk */
405        if (rank == 0) {
406          strcpy(dataNameStr, "escript_scalar_data");
407          sprintf(dataTypeStr, "Scalars=\"%s\"", dataNameStr);
408          nCompReqd = 1;
409        }
410        else if (rank == 1) {
411          strcpy(dataNameStr, "escript_vector_data");
412          sprintf(dataTypeStr, "Vectors=\"%s\"", dataNameStr);
413          nCompReqd = 3;
414        }
415        else if (rank == 2) {
416          strcpy(dataNameStr, "escript_tensor_data");
417          sprintf(dataTypeStr, "Tensors=\"%s\"", dataNameStr);
418          nCompReqd = 9;
419        }
420        /* if have cell centred data then use CellData element,
421         * if have node centred data, then use PointData element
422         */
423        if (isCellCentered) {
424          /* now for the cell data */
425          fprintf(fileHandle_p, "<CellData %s>\n", dataTypeStr);
426          fprintf(fileHandle_p,
427              "<DataArray "
428              "Name=\"%s\" "
429              "type=\"Float32\" "
430              "NumberOfComponents=\"%d\" "
431              "format=\"ascii\">\n",
432              dataNameStr, nCompReqd);
433          int numPointsPerSample = elements->ReferenceElement->numQuadNodes;
434          if (numPointsPerSample) {
435        for (i=0; i<numCells; i++) {
436          values = getSampleData(data_p, i);
437          double sampleAvg[nComp];
438          for (k=0; k<nComp; k++) {
439            /* averaging over the number of points in the sample */
440            rtmp = 0.;
441            for (j=0; j<numPointsPerSample; j++) {
442              rtmp += values[INDEX2(k,j,nComp)];
443            }
444            sampleAvg[k] = rtmp/numPointsPerSample;
445          }
446          /* if the number of required components is more than the number
447           * of actual components, pad with zeros
448           */
449          /* probably only need to get shape of first element */
450          int shape = getDataPointShape(data_p, 0);
451          if (shape > 3) {
452            Finley_ErrorCode = VALUE_ERROR;
453            sprintf(Finley_ErrorMsg,
454                "shape should be 1, 2, or 3; I got %d\n", shape);
455            return;
456          }
457          /* write the data different ways for scalar, vector and tensor */
458          if (nCompReqd == 1) {
459            fprintf(fileHandle_p, " %f", sampleAvg[0]);
460          }
461          else if (nCompReqd == 3) {
462            int shape = getDataPointShape(data_p, 0);
463            /* write out the data */
464            for (int m=0; m<shape; m++) {
465              fprintf(fileHandle_p, " %f", sampleAvg[m]);
466            }
467            /* pad with zeros */
468            for (int m=0; m<nCompReqd-shape; m++) {
469              fprintf(fileHandle_p, " 0");
470            }
471          }
472          else if (nCompReqd == 9) {
473            /* tensor data, so have a 3x3 matrix to output as a row
474             * of 9 data points */
475            int count = 0;
476            int elems = 0;
477            for (int m=0; m<shape; m++) {
478              for (int n=0; n<shape; n++) {
479            fprintf(fileHandle_p, " %f", sampleAvg[count]);
480            count++;
481            elems++;
482              }
483              for (int n=0; n<3-shape; n++) {
484            fprintf(fileHandle_p, " 0");
485            elems++;
486              }
487            }
488            for (int m=0; m<nCompReqd-elems; m++) {
489              fprintf(fileHandle_p, " 0");
490            }
491          }
492          fprintf(fileHandle_p, "\n");
493        }
494        }        }
495        if (isCellCentered) {        fprintf(fileHandle_p, "</DataArray>\n");
496            int numPointsPerSample=elements->ReferenceElement->numQuadNodes;        fprintf(fileHandle_p, "</CellData>\n");
497            if (numPointsPerSample) {      } else {
498               fprintf(fileHandle_p, "items %d data follows\n", elements->numElements);        /* now for the point data */
499               for (i=0;i<elements->numElements;i++) {        fprintf(fileHandle_p, "<PointData %s>\n", dataTypeStr);
500                   values=getSampleData(data_p,i);        fprintf(fileHandle_p, "<DataArray "
501                   for (k=0;k<nComp;k++) {            "Name=\"%s\" "
502                       rtmp=0.;            "type=\"Float32\" "
503                       for (j=0;j<numPointsPerSample;j++) rtmp+=values[INDEX2(k,j,nComp)];            "NumberOfComponents=\"%d\" "
504                       fprintf(fileHandle_p, " %f", rtmp/numPointsPerSample);            "format=\"ascii\">\n",
505                   }            dataNameStr, nComp);
506               fprintf(fileHandle_p, "\n");        for (i=0; i<mesh_p->Nodes->numNodes; i++) {
507               }      switch (nodetype) {
508               fprintf(fileHandle_p, "attribute \"dep\" string \"connections\"\n");      case(FINLEY_DEGREES_OF_FREEDOM):
509           }        values = getSampleData(data_p,
510        } else {                   mesh_p->Nodes->degreeOfFreedom[i]);
511            fprintf(fileHandle_p, "items %d data follows\n", mesh_p->Nodes->reducedNumNodes);        break;
512            for (i=0;i<mesh_p->Nodes->numNodes;i++) {      case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):
513                if (mesh_p->Nodes->toReduced[i]>=0) {        if (mesh_p->Nodes->toReduced[i]>=0) {
514                   switch (nodetype) {          values = getSampleData(data_p,
515                      case(FINLEY_DEGREES_OF_FREEDOM):                     mesh_p->Nodes->reducedDegreeOfFreedom[i]);
516                         values=getSampleData(data_p,mesh_p->Nodes->degreeOfFreedom[i]);        }
517                         break;        break;
518                      case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):      case(FINLEY_NODES):
519                         values=getSampleData(data_p,mesh_p->Nodes->reducedDegreeOfFreedom[i]);        values = getSampleData(data_p,i);
520                         break;        break;
521                      case(FINLEY_NODES):      }
522                         values=getSampleData(data_p,i);      /* write out the data */
523                         break;      /* if the number of required components is more than the number
524                   }       * of actual components, pad with zeros
525                }       */
526                for (k=0;k<nComp;k++) fprintf(fileHandle_p, " %f", values[k]);      /* probably only need to get shape of first element */
527            fprintf(fileHandle_p, "\n");      int shape = getDataPointShape(data_p, 0);
528            }      if (shape > 3) {
529            fprintf(fileHandle_p, "attribute \"dep\" string \"positions\"\n");        Finley_ErrorCode = VALUE_ERROR;
530          sprintf(Finley_ErrorMsg,
531              "shape should be 1, 2, or 3; I got %d\n", shape);
532          return;
533        }
534        /* write the data different ways for scalar, vector and tensor */
535        if (nCompReqd == 1) {
536          fprintf(fileHandle_p, " %f", values[0]);
537        }
538        else if (nCompReqd == 3) {
539          int shape = getDataPointShape(data_p, 0);
540          /* write out the data */
541          for (int m=0; m<shape; m++) {
542            fprintf(fileHandle_p, " %f", values[m]);
543          }
544          /* pad with zeros */
545          for (int m=0; m<nCompReqd-shape; m++) {
546            fprintf(fileHandle_p, " 0");
547          }
548        }
549        else if (nCompReqd == 9) {
550          /* tensor data, so have a 3x3 matrix to output as a row
551           * of 9 data points */
552          int count = 0;
553          int elems = 0;
554          for (int m=0; m<shape; m++) {
555            for (int n=0; n<shape; n++) {
556              fprintf(fileHandle_p, " %f", values[count]);
557              count++;
558              elems++;
559            }
560            for (int n=0; n<3-shape; n++) {
561              fprintf(fileHandle_p, " 0");
562              elems++;
563            }
564          }
565          for (int m=0; m<nCompReqd-elems; m++) {
566            fprintf(fileHandle_p, " 0");
567          }
568        }
569        fprintf(fileHandle_p, "\n");
570        }        }
571          fprintf(fileHandle_p, "</DataArray>\n");
572          fprintf(fileHandle_p, "</PointData>\n");
573        }
574    }    }
575    
576    /* and finish it up */    /* finish off the piece */
577    fprintf(fileHandle_p, "object 4 class field\n");    fprintf(fileHandle_p, "</Piece>\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");  
578    
579    fprintf(fileHandle_p, "</UnstructuredGrid>\n");    fprintf(fileHandle_p, "</UnstructuredGrid>\n");
580    /* write the xml footer */    /* write the xml footer */
# Line 234  void Finley_Mesh_saveVTK(const char * fi Line 583  void Finley_Mesh_saveVTK(const char * fi
583    fclose(fileHandle_p);    fclose(fileHandle_p);
584    return;    return;
585  }  }
   
 /*  
  * $Log$  
  * Revision 1.2  2005/02/14 04:14:42  jgs  
  * *** empty log message ***  
  *  
  * Revision 1.1.2.2  2005/02/10 01:34:22  cochrane  
  * Quick fix to make sure that saveVTK compiles so that finley is still buildable.  Apologies to those this has affected.  
  *  
  * Revision 1.1.2.1  2005/02/09 06:53:15  cochrane  
  * Initial import to repository.  This is the file to implement saving finley/escript meshes out to vtk formatted files.  It is basically just a hack of the opendx equivalent, with a lot of the opendx stuff still in the file, so it doesn't actually work just yet, but it probably needs to be added to the cvs.  
  *  
  * Revision 1.1.1.1  2004/10/26 06:53:57  jgs  
  * initial import of project esys2  
  *  
  * Revision 1.1  2004/07/27 08:27:11  gross  
  * Finley: saveDX added: now it is possible to write data on boundary and contact elements  
  *  
  */  

Legend:
Removed from v.112  
changed lines
  Added in v.113

  ViewVC Help
Powered by ViewVC 1.1.26