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

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

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

revision 110 by jgs, Mon Feb 14 04:14:42 2005 UTC revision 121 by jgs, Fri May 6 04:26:16 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        int counter = 0;
352        for (i = 0; i < numCells; i++) {
353          fprintf(fileHandle_p, "%d ",
354              mesh_p->Elements->Nodes[INDEX2(0, i, NN)]);
355          counter++; /* counter for the number of connectivity points written */
356          /* if the counter gets too big (i.e. the line gets too long),
357           * then add a newline and reset */
358          if (counter > 19) {
359          fprintf(fileHandle_p, "\n");
360          counter = 0;
361          }
362          for (j = 1; j < numVTKNodesPerElement; j++) {
363        fprintf(fileHandle_p,"%d ", mesh_p->Elements->Nodes[INDEX2(j, i, NN)]);
364        counter++;
365        /* if the counter gets too big (i.e. the line gets too long),
366         * then add a newline and reset */
367        if (counter > 19) {
368            fprintf(fileHandle_p, "\n");
369            counter = 0;
370        }
371          }
372        }
373        fprintf(fileHandle_p, "\n");
374        fprintf(fileHandle_p, "</DataArray>\n");
375    
376        /* write out the DataArray element for the offsets */
377        fprintf(fileHandle_p, "<DataArray "
378            "Name=\"offsets\" "
379            "type=\"Int32\" "
380            "format=\"ascii\">\n");
381        counter = 0;  /* counter for the number of offsets written to file */
382        for (i=numVTKNodesPerElement; i<=numCells*numVTKNodesPerElement; i+=numVTKNodesPerElement) {
383          fprintf(fileHandle_p, "%d ", i);
384          counter++;
385          /* if the counter gets too big (i.e. the line gets too long),
386           * then add a newline and reset */
387          if (counter > 19) {
388          counter = 0;
389          fprintf(fileHandle_p, "\n");
390          }
391        }
392        fprintf(fileHandle_p, "\n");
393        fprintf(fileHandle_p, "</DataArray>\n");
394    
395    /* positions */      /* write out the DataArray element for the types */
396    fprintf(fileHandle_p, "object 1 class array type float rank 1 shape %d items %d data follows\n", nDim, mesh_p->Nodes->reducedNumNodes);      counter = 0; /* counter for the number of types written to file */
397    for (i = 0; i < mesh_p->Nodes->numNodes; i++) {      fprintf(fileHandle_p, "<DataArray "
398      if (mesh_p->Nodes->toReduced[i]>=0) {          "Name=\"types\" "
399         fprintf(fileHandle_p, "%e", mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)]);          "type=\"UInt8\" "
400         for (j = 1; j < nDim; j++) fprintf(fileHandle_p, " %f",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);          "format=\"ascii\">\n");
401         fprintf(fileHandle_p, "\n");      for (i=0; i<numCells; i++) {
402          fprintf(fileHandle_p, "%d ", cellType);
403          counter++;
404          /* if the counter gets too big (i.e. the line gets too long),
405           * then add a newline and reset */
406          if (counter > 30) {
407          counter = 0;
408          fprintf(fileHandle_p, "\n");
409          }
410      }      }
411    }      fprintf(fileHandle_p, "\n");
412    /* 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");  
413    
414        /* finish off the <Cells> element */
415        fprintf(fileHandle_p, "</Cells>\n");
416    }    }
417    
418    /* data */    /* data */
419    if (!isEmpty(data_p)) {    if (!isEmpty(data_p)) {
420        int rank=getDataPointRank(data_p);      int rank = getDataPointRank(data_p);
421        int* shape=getDataPointShape(data_p);      int nComp = getDataPointSize(data_p);
422        int nComp=getDataPointSize(data_p);      /* barf if rank is greater than two */
423        fprintf(fileHandle_p, "object 3 class array type float rank %d ", rank);      if (rank > 2) {
424        if (0 < rank) {        Finley_ErrorCode = VALUE_ERROR;
425           fprintf(fileHandle_p, "shape ");        sprintf(Finley_ErrorMsg,
426           for (i = 0; i < rank; i++) fprintf(fileHandle_p, "%d ", shape[i]);            "Vtk can't handle objects with rank greater than 2\n"
427        }            "object rank = %d\n", rank);
428        if (isCellCentered) {        return;
429            int numPointsPerSample=elements->ReferenceElement->numQuadNodes;      }
430            if (numPointsPerSample) {      /* if the rank == 0:   --> scalar data
431               fprintf(fileHandle_p, "items %d data follows\n", elements->numElements);       * if the rank == 1:   --> vector data
432               for (i=0;i<elements->numElements;i++) {       * if the rank == 2:   --> tensor data
433                   values=getSampleData(data_p,i);       */
434                   for (k=0;k<nComp;k++) {      char dataNameStr[31], dataTypeStr[63];
435                       rtmp=0.;      int nCompReqd;   /* the number of components required by vtk */
436                       for (j=0;j<numPointsPerSample;j++) rtmp+=values[INDEX2(k,j,nComp)];      if (rank == 0) {
437                       fprintf(fileHandle_p, " %f", rtmp/numPointsPerSample);        strcpy(dataNameStr, "escript_scalar_data");
438                   }        sprintf(dataTypeStr, "Scalars=\"%s\"", dataNameStr);
439               fprintf(fileHandle_p, "\n");        nCompReqd = 1;
440               }      }
441               fprintf(fileHandle_p, "attribute \"dep\" string \"connections\"\n");      else if (rank == 1) {
442           }        strcpy(dataNameStr, "escript_vector_data");
443        } else {        sprintf(dataTypeStr, "Vectors=\"%s\"", dataNameStr);
444            fprintf(fileHandle_p, "items %d data follows\n", mesh_p->Nodes->reducedNumNodes);        nCompReqd = 3;
445            for (i=0;i<mesh_p->Nodes->numNodes;i++) {      }
446                if (mesh_p->Nodes->toReduced[i]>=0) {      else if (rank == 2) {
447                   switch (nodetype) {        strcpy(dataNameStr, "escript_tensor_data");
448                      case(FINLEY_DEGREES_OF_FREEDOM):        sprintf(dataTypeStr, "Tensors=\"%s\"", dataNameStr);
449                         values=getSampleData(data_p,mesh_p->Nodes->degreeOfFreedom[i]);        nCompReqd = 9;
450                         break;      }
451                      case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):      /* if have cell centred data then use CellData element,
452                         values=getSampleData(data_p,mesh_p->Nodes->reducedDegreeOfFreedom[i]);       * if have node centred data, then use PointData element
453                         break;       */
454                      case(FINLEY_NODES):      if (isCellCentered) {
455                         values=getSampleData(data_p,i);        /* now for the cell data */
456                         break;        fprintf(fileHandle_p, "<CellData %s>\n", dataTypeStr);
457                   }        fprintf(fileHandle_p,
458                }            "<DataArray "
459                for (k=0;k<nComp;k++) fprintf(fileHandle_p, " %f", values[k]);            "Name=\"%s\" "
460            fprintf(fileHandle_p, "\n");            "type=\"Float32\" "
461            }            "NumberOfComponents=\"%d\" "
462            fprintf(fileHandle_p, "attribute \"dep\" string \"positions\"\n");            "format=\"ascii\">\n",
463              dataNameStr, nCompReqd);
464          int numPointsPerSample = elements->ReferenceElement->numQuadNodes;
465          if (numPointsPerSample) {
466        for (i=0; i<numCells; i++) {
467          values = getSampleData(data_p, i);
468          double sampleAvg[nComp];
469          for (k=0; k<nComp; k++) {
470            /* averaging over the number of points in the sample */
471            rtmp = 0.;
472            for (j=0; j<numPointsPerSample; j++) {
473              rtmp += values[INDEX2(k,j,nComp)];
474            }
475            sampleAvg[k] = rtmp/numPointsPerSample;
476          }
477          /* if the number of required components is more than the number
478           * of actual components, pad with zeros
479           */
480          /* probably only need to get shape of first element */
481          int shape = getDataPointShape(data_p, 0);
482          if (shape > 3) {
483            Finley_ErrorCode = VALUE_ERROR;
484            sprintf(Finley_ErrorMsg,
485                "shape should be 1, 2, or 3; I got %d\n", shape);
486            return;
487          }
488          /* write the data different ways for scalar, vector and tensor */
489          if (nCompReqd == 1) {
490            fprintf(fileHandle_p, " %f", sampleAvg[0]);
491          }
492          else if (nCompReqd == 3) {
493            int shape = getDataPointShape(data_p, 0);
494            /* write out the data */
495            for (int m=0; m<shape; m++) {
496              fprintf(fileHandle_p, " %f", sampleAvg[m]);
497            }
498            /* pad with zeros */
499            for (int m=0; m<nCompReqd-shape; m++) {
500              fprintf(fileHandle_p, " 0");
501            }
502          }
503          else if (nCompReqd == 9) {
504            /* tensor data, so have a 3x3 matrix to output as a row
505             * of 9 data points */
506            int count = 0;
507            int elems = 0;
508            for (int m=0; m<shape; m++) {
509              for (int n=0; n<shape; n++) {
510            fprintf(fileHandle_p, " %f", sampleAvg[count]);
511            count++;
512            elems++;
513              }
514              for (int n=0; n<3-shape; n++) {
515            fprintf(fileHandle_p, " 0");
516            elems++;
517              }
518            }
519            for (int m=0; m<nCompReqd-elems; m++) {
520              fprintf(fileHandle_p, " 0");
521            }
522          }
523          fprintf(fileHandle_p, "\n");
524        }
525          }
526          fprintf(fileHandle_p, "</DataArray>\n");
527          fprintf(fileHandle_p, "</CellData>\n");
528        } else {
529          /* now for the point data */
530          fprintf(fileHandle_p, "<PointData %s>\n", dataTypeStr);
531          fprintf(fileHandle_p, "<DataArray "
532              "Name=\"%s\" "
533              "type=\"Float32\" "
534              "NumberOfComponents=\"%d\" "
535              "format=\"ascii\">\n",
536              dataNameStr, nCompReqd);
537          for (i=0; i<mesh_p->Nodes->numNodes; i++) {
538        switch (nodetype) {
539        case(FINLEY_DEGREES_OF_FREEDOM):
540          values = getSampleData(data_p,
541                     mesh_p->Nodes->degreeOfFreedom[i]);
542          break;
543        case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):
544          if (mesh_p->Nodes->toReduced[i]>=0) {
545            values = getSampleData(data_p,
546                       mesh_p->Nodes->reducedDegreeOfFreedom[i]);
547          }
548          break;
549        case(FINLEY_NODES):
550          values = getSampleData(data_p,i);
551          break;
552        }
553        /* write out the data */
554        /* if the number of required components is more than the number
555         * of actual components, pad with zeros
556         */
557        /* probably only need to get shape of first element */
558        int shape = getDataPointShape(data_p, 0);
559        if (shape > 3) {
560          Finley_ErrorCode = VALUE_ERROR;
561          sprintf(Finley_ErrorMsg,
562              "shape should be 1, 2, or 3; I got %d\n", shape);
563          return;
564        }
565        /* write the data different ways for scalar, vector and tensor */
566        if (nCompReqd == 1) {
567          fprintf(fileHandle_p, " %f", values[0]);
568        }
569        else if (nCompReqd == 3) {
570          int shape = getDataPointShape(data_p, 0);
571          /* write out the data */
572          for (int m=0; m<shape; m++) {
573            fprintf(fileHandle_p, " %f", values[m]);
574          }
575          /* pad with zeros */
576          for (int m=0; m<nCompReqd-shape; m++) {
577            fprintf(fileHandle_p, " 0");
578          }
579        }
580        else if (nCompReqd == 9) {
581          /* tensor data, so have a 3x3 matrix to output as a row
582           * of 9 data points */
583          int count = 0;
584          int elems = 0;
585          for (int m=0; m<shape; m++) {
586            for (int n=0; n<shape; n++) {
587              fprintf(fileHandle_p, " %f", values[count]);
588              count++;
589              elems++;
590            }
591            for (int n=0; n<3-shape; n++) {
592              fprintf(fileHandle_p, " 0");
593              elems++;
594            }
595          }
596          for (int m=0; m<nCompReqd-elems; m++) {
597            fprintf(fileHandle_p, " 0");
598          }
599        }
600        fprintf(fileHandle_p, "\n");
601        }        }
602          fprintf(fileHandle_p, "</DataArray>\n");
603          fprintf(fileHandle_p, "</PointData>\n");
604        }
605    }    }
606    
607    /* and finish it up */    /* finish off the piece */
608    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");  
609    
610    fprintf(fileHandle_p, "</UnstructuredGrid>\n");    fprintf(fileHandle_p, "</UnstructuredGrid>\n");
611    /* write the xml footer */    /* write the xml footer */
# Line 237  void Finley_Mesh_saveVTK(const char * fi Line 617  void Finley_Mesh_saveVTK(const char * fi
617    
618  /*  /*
619   * $Log$   * $Log$
620   * Revision 1.2  2005/02/14 04:14:42  jgs   * Revision 1.4  2005/05/06 04:26:15  jgs
621   * *** empty log message ***   * Merge of development branch back to main trunk on 2005-05-06
622     *
623     * Revision 1.1.2.6  2005/05/06 01:17:19  cochrane
624     * Fixed incorrect reporting of number of components in PointData arrays for
625     * vector data.
626     *
627     * Revision 1.1.2.5  2005/05/05 05:38:44  cochrane
628     * Improved formatting of VTK file output.
629     *
630     * Revision 1.1.2.4  2005/02/22 10:03:54  cochrane
631     * Implementation of writing of vtk xml files from finley.  This function will
632     * require more testing, but on the cases that I have tried (and with the help
633     * of Lutz and mayavi), it looks like it's producing the correct output.  Testing
634     * with more realistic data would be good.  I'm at least confident enough to
635     * commit my changes.
636     *
637     * Revision 1.1.2.3  2005/02/17 05:53:26  gross
638     * some bug in saveDX fixed: in fact the bug was in
639     * DataC/getDataPointShape
640   *   *
641   * Revision 1.1.2.2  2005/02/10 01:34:22  cochrane   * Revision 1.1.2.2  2005/02/10 01:34:22  cochrane
642   * Quick fix to make sure that saveVTK compiles so that finley is still buildable.  Apologies to those this has affected.   * Quick fix to make sure that saveVTK compiles so that finley is still buildable.  Apologies to those this has affected.

Legend:
Removed from v.110  
changed lines
  Added in v.121

  ViewVC Help
Powered by ViewVC 1.1.26