/[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 121 by jgs, Fri May 6 04:26:16 2005 UTC revision 153 by jgs, Tue Oct 25 01:51:20 2005 UTC
# Line 1  Line 1 
1  /* $Id$ */  /*
2     ******************************************************************************
3     *                                                                            *
4     *       COPYRIGHT  ACcESS 2003,2004,2005 -  All Rights Reserved              *
5     *                                                                            *
6     * This software is the property of ACcESS. No part of this code              *
7     * may be copied in any form or by any means without the expressed written    *
8     * consent of ACcESS.  Copying, use or modification of this software          *
9     * by any unauthorised person is illegal unless that person has a software    *
10     * license agreement with ACcESS.                                             *
11     *                                                                            *
12     ******************************************************************************
13    */
14    
15    
16  /**************************************************************/  /**************************************************************/
17    
# Line 6  Line 19 
19    
20  /**************************************************************/  /**************************************************************/
21    
 /*   Copyrights by ACcESS, Australia 2004 */  
22  /*   Author: Paul Cochrane, cochrane@esscc.uq.edu.au */  /*   Author: Paul Cochrane, cochrane@esscc.uq.edu.au */
23    /*   Version: $Id$ */
24    
25  /**************************************************************/  /**************************************************************/
26    
 #include "Finley.h"  
 #include "Common.h"  
27  #include "Mesh.h"  #include "Mesh.h"
 #include "escript/Data/DataC.h"  
28  #include "vtkCellType.h"  /* copied from vtk source directory !!! */  #include "vtkCellType.h"  /* copied from vtk source directory !!! */
29    
30  void Finley_Mesh_saveVTK(const char * filename_p, Finley_Mesh *mesh_p, escriptDataC* data_p) {  /**************************************************************/
31    
32    void Finley_Mesh_saveVTK(const char * filename_p, Finley_Mesh *mesh_p, const dim_t num_data,char* *names_p, escriptDataC* *data_pp) {
33      char error_msg[LenErrorMsg_MAX];
34    /* if there is no mesh we just return */    /* if there is no mesh we just return */
35    if (mesh_p==NULL) return;    if (mesh_p==NULL) return;
36    Finley_ElementFile* elements=NULL;  
37    char elemTypeStr[32];    int i, j, k, numVTKNodesPerElement,i_data;
38    int i, j, k, numVTKNodesPerElement, isCellCentered, nodetype;    index_t j2;
39    double* values, rtmp;    double* values, rtmp;
40    int nDim = mesh_p->Nodes->numDim;    
41      /* open the file and check handle */
42    
43    /* get a pointer to the relevant mesh component */    FILE * fileHandle_p = fopen(filename_p, "w");
44    if (isEmpty(data_p)) {    if (fileHandle_p==NULL) {
45      elements=mesh_p->Elements;      sprintf(error_msg, "saveVTK: File %s could not be opened for writing.", filename_p);
46        Finley_setError(IO_ERROR,error_msg);
47        return;
48      }
49      /* find the mesh type to be written */
50      int nodetype=FINLEY_DEGREES_OF_FREEDOM;
51      int elementtype=FINLEY_UNKNOWN;
52      bool_t isCellCentered[num_data],write_celldata=FALSE,write_pointdata=FALSE;
53      for (i_data=0;i_data<num_data;++i_data) {
54         if (! isEmpty(data_pp[i_data])) {
55            switch(getFunctionSpaceType(data_pp[i_data])) {
56               case FINLEY_DEGREES_OF_FREEDOM:
57                 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
58                 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
59                     elementtype=FINLEY_ELEMENTS;
60                 } else {
61                     Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
62                     fclose(fileHandle_p);
63                     return;
64                 }
65                 isCellCentered[i_data]=FALSE;
66                 break;
67               case FINLEY_REDUCED_DEGREES_OF_FREEDOM:
68                 nodetype = FINLEY_REDUCED_DEGREES_OF_FREEDOM;
69                 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
70                     elementtype=FINLEY_ELEMENTS;
71                 } else {
72                     Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
73                     fclose(fileHandle_p);
74                     return;
75                 }
76                 isCellCentered[i_data]=FALSE;
77                 break;
78               case FINLEY_NODES:
79                 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
80                 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
81                     elementtype=FINLEY_ELEMENTS;
82                 } else {
83                     Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
84                     fclose(fileHandle_p);
85                     return;
86                 }
87                 isCellCentered[i_data]=FALSE;
88                 break;
89               case FINLEY_ELEMENTS:
90                 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
91                 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
92                     elementtype=FINLEY_ELEMENTS;
93                 } else {
94                     Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
95                     fclose(fileHandle_p);
96                     return;
97                 }
98                 isCellCentered[i_data]=TRUE;
99                 break;
100               case FINLEY_FACE_ELEMENTS:
101                 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
102                 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_FACE_ELEMENTS) {
103                     elementtype=FINLEY_FACE_ELEMENTS;
104                 } else {
105                     Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
106                     fclose(fileHandle_p);
107                     return;
108                 }
109                 isCellCentered[i_data]=TRUE;
110                 break;
111               case FINLEY_POINTS:
112                 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
113                 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_POINTS) {
114                     elementtype=FINLEY_POINTS;
115                 } else {
116                     Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
117                     fclose(fileHandle_p);
118                     return;
119                 }
120                 isCellCentered[i_data]=TRUE;
121                 break;
122               case FINLEY_CONTACT_ELEMENTS_1:
123                 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
124                 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {
125                     elementtype=FINLEY_CONTACT_ELEMENTS_1;
126                 } else {
127                     Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
128                     fclose(fileHandle_p);
129                     return;
130                 }
131                 isCellCentered[i_data]=TRUE;
132                 break;
133               case FINLEY_CONTACT_ELEMENTS_2:
134                 nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
135                 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {
136                     elementtype=FINLEY_CONTACT_ELEMENTS_1;
137                 } else {
138                     Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
139                     fclose(fileHandle_p);
140                     return;
141                 }
142                 isCellCentered[i_data]=TRUE;
143                 break;
144               default:
145                 sprintf(error_msg,"saveVTK: Finley does not know anything about function space type %d",getFunctionSpaceType(data_pp[i_data]));
146                 Finley_setError(TYPE_ERROR,error_msg);
147                 fclose(fileHandle_p);
148                 return;
149            }
150            if (isCellCentered[i_data]) {
151               write_celldata=TRUE;
152            } else {
153               write_pointdata=TRUE;
154            }
155         }
156      }
157      /* select nomber of points and the mesh component */
158      int numPoints = mesh_p->Nodes->numNodes;
159      if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
160          numPoints = mesh_p->Nodes->reducedNumNodes;
161    } else {    } else {
162      switch(getFunctionSpaceType(data_p)) {        numPoints = mesh_p->Nodes->numNodes;
163      case(FINLEY_DEGREES_OF_FREEDOM):    }
164        nodetype = FINLEY_DEGREES_OF_FREEDOM;    if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;
165        isCellCentered = FALSE;    Finley_ElementFile* elements=NULL;
166        elements = mesh_p->Elements;    switch(elementtype) {
167        break;      case FINLEY_ELEMENTS:
     case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):  
       Finley_ErrorCode=VALUE_ERROR;  
       sprintf(Finley_ErrorMsg,  
           "Reduced degrees of freedom is not yet "  
           "implemented for saving vtk files\n");  
       return;  
     case(FINLEY_NODES):  
       nodetype=FINLEY_NODES;  
       isCellCentered=FALSE;  
       elements=mesh_p->Elements;  
       break;  
     case(FINLEY_ELEMENTS):  
       isCellCentered=TRUE;  
168        elements=mesh_p->Elements;        elements=mesh_p->Elements;
169        break;        break;
170      case(FINLEY_FACE_ELEMENTS):      case FINLEY_FACE_ELEMENTS:
       isCellCentered=TRUE;  
171        elements=mesh_p->FaceElements;        elements=mesh_p->FaceElements;
172        break;        break;
173      case(FINLEY_POINTS):      case FINLEY_POINTS:
       isCellCentered=TRUE;  
174        elements=mesh_p->Points;        elements=mesh_p->Points;
175        break;        break;
176      case(FINLEY_CONTACT_ELEMENTS_1):      case FINLEY_CONTACT_ELEMENTS_1:
     case(FINLEY_CONTACT_ELEMENTS_2):  
       isCellCentered=TRUE;  
177        elements=mesh_p->ContactElements;        elements=mesh_p->ContactElements;
178        break;        break;
     default:  
       Finley_ErrorCode=TYPE_ERROR;  
       sprintf(Finley_ErrorMsg,  
           "Finley does not know anything about function space type %d",  
           getFunctionSpaceType(data_p));  
       return;  
     }  
179    }    }
180      if (elements==NULL) {
181    /* the number of points */       Finley_setError(SYSTEM_ERROR,"saveVTK: undefined element file");
182    int numPoints = mesh_p->Nodes->numNodes;       fclose(fileHandle_p);
183         return;
   /* the number of cells */  
   if (elements == NULL) {  
     Finley_ErrorCode = VALUE_ERROR;  
     sprintf(Finley_ErrorMsg,  
         "elements object is NULL; cannot proceed");  
     return;  
184    }    }
185      /* map finley element type to VTK element type */
186    int numCells = elements->numElements;      int numCells = elements->numElements;  
187        int cellType;
188    /* open the file and check handle */    ElementTypeId TypeId;
189    FILE * fileHandle_p = fopen(filename_p, "w");    char elemTypeStr[32];
190    if (fileHandle_p==NULL) {    if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
191      Finley_ErrorCode=IO_ERROR;       TypeId = elements->LinearReferenceElement->Type->TypeId;
     sprintf(Finley_ErrorMsg,  
         "File %s could not be opened for writing.", filename_p);  
     return;  
   }  
   /* xml header */  
   fprintf(fileHandle_p, "<?xml version=\"1.0\"?>\n");  
   fprintf(fileHandle_p,  
       "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\">\n");  
   
   /* finley uses an unstructured mesh, so UnstructuredGrid *should* work */  
   fprintf(fileHandle_p, "<UnstructuredGrid>\n");  
   
   /* is there only one "piece" to the data?? */  
   fprintf(fileHandle_p, "<Piece "  
       "NumberOfPoints=\"%d\" "  
       "NumberOfCells=\"%d\">\n",  
       numPoints, numCells);  
   
   /* now for the points; equivalent to positions section in saveDX() */  
   /* "The points element explicitly defines coordinates for each point  
    * individually.  It contains one DataArray element describing an array  
    * with three components per value, each specifying the coordinates of one  
    * point" - from Vtk User's Guide  
    */  
   fprintf(fileHandle_p, "<Points>\n");  
   /*  
    * the reason for this if statement is explained in the long comment below  
    */  
   if (nDim < 3) {  
     fprintf(fileHandle_p, "<DataArray "  
         "NumberOfComponents=\"3\" "  
         "type=\"Float32\" "  
         "format=\"ascii\">\n");  
192    } else {    } else {
193      fprintf(fileHandle_p, "<DataArray "       TypeId = elements->ReferenceElement->Type->TypeId;
         "NumberOfComponents=\"%d\" "  
         "type=\"Float32\" "  
         "format=\"ascii\">\n",  
         nDim);  
194    }    }
   for (i = 0; i < mesh_p->Nodes->numNodes; i++) {  
     fprintf(fileHandle_p,  
         "%e", mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)]);  
     for (j = 1; j < nDim; j++) {  
       fprintf(fileHandle_p,  
           " %f",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);  
       /* vtk/mayavi doesn't like 2D data, it likes 3D data with a degenerate  
        * third dimension to handle 2D data (like a sheet of paper).  So, if  
        * nDim is 2, we have to append zeros to the array to get this third  
        * dimension, and keep the visualisers happy.  
        * Indeed, if nDim is less than 3, must pad all empty dimensions, so  
        * that the total number of dims is 3.  
        */  
       if (nDim < 3) {  
     for (k=0; k<3-nDim; k++) {  
       fprintf(fileHandle_p, " 0");  
     }  
       }  
     }  
     fprintf(fileHandle_p, "\n");  
   }  
   fprintf(fileHandle_p, "</DataArray>\n");  
   fprintf(fileHandle_p, "</Points>\n");  
195    
196    /* connections */    switch(TypeId) {
   /* now for the cells */  
   /* "The Cells element defines cells explicitly by specifying point  
    * connectivity and cell types.  It contains three DataArray elements.  The  
    * first array specifies the point connectivity.  All cells' point lists  
    * are concatenated together.  The second array specifies th eoffset into  
    * the connectivity array for the end of each cell.  The third array  
    * specifies the type of each cell.  
    */  
   /* if no element table is present jump over the connection table */  
   int cellType;  
   if (elements!=NULL) {  
     fprintf(fileHandle_p, "<Cells>\n");  
     ElementTypeId TypeId = elements->ReferenceElement->Type->TypeId;  
     switch(TypeId) {  
197      case Point1:      case Point1:
       cellType = VTK_VERTEX;  
       break;  
     case Line2:  
       cellType = VTK_LINE;  
       break;  
     case Line3:  
       cellType = VTK_QUADRATIC_EDGE;  
       break;  
     case Tri3:  
       cellType = VTK_TRIANGLE;  
       break;  
     case Tri6:  
       cellType = VTK_QUADRATIC_TRIANGLE;  
       break;  
     case Rec4:  
       cellType = VTK_QUAD;  
       break;  
     case Rec8:  
       cellType = VTK_QUADRATIC_QUAD;  
       break;  
     case Tet4:  
       cellType = VTK_TETRA;  
       break;  
     case Tet10:  
       cellType = VTK_QUADRATIC_TETRA;  
       break;  
     case Hex8:  
       cellType = VTK_HEXAHEDRON;  
       break;  
     case Hex20:  
       cellType = VTK_QUADRATIC_HEXAHEDRON;  
       break;  
198      case Line2Face:      case Line2Face:
       cellType = VTK_VERTEX;  
       break;  
199      case Line3Face:      case Line3Face:
200        case Point1_Contact:
201        case Line2Face_Contact:
202        case Line3Face_Contact:
203        cellType = VTK_VERTEX;        cellType = VTK_VERTEX;
204          numVTKNodesPerElement = 1;
205          strcpy(elemTypeStr, "VTK_VERTEX");
206        break;        break;
207    
208        case Line2:
209      case Tri3Face:      case Tri3Face:
       cellType = VTK_LINE;  
       break;  
     case Tri6Face:  
       cellType = VTK_QUADRATIC_EDGE;  
       break;  
210      case Rec4Face:      case Rec4Face:
       cellType = VTK_LINE;  
       break;  
     case Rec8Face:  
       cellType = VTK_QUADRATIC_EDGE;  
       break;  
     case Tet4Face:  
       cellType = VTK_TRIANGLE;  
       break;  
     case Tet10Face:  
       cellType = VTK_QUADRATIC_TRIANGLE;  
       break;  
     case Hex8Face:  
       cellType = VTK_QUADRATIC_QUAD;  
       break;  
     case Hex20Face:  
       cellType = VTK_QUADRATIC_QUAD;  
       break;  
     case Point1_Contact:  
       cellType = VTK_VERTEX;  
       break;  
211      case Line2_Contact:      case Line2_Contact:
       cellType = VTK_LINE;  
       break;  
     case Line3_Contact:  
       cellType = VTK_QUADRATIC_EDGE;  
       break;  
212      case Tri3_Contact:      case Tri3_Contact:
       cellType = VTK_TRIANGLE;  
       break;  
     case Tri6_Contact:  
       cellType = VTK_QUADRATIC_TRIANGLE;  
       break;  
     case Rec4_Contact:  
       cellType = VTK_QUAD;  
       break;  
     case Rec8_Contact:  
       cellType = VTK_QUADRATIC_QUAD;  
       break;  
     case Line2Face_Contact:  
       cellType = VTK_VERTEX;  
       break;  
     case Line3Face_Contact:  
       cellType = VTK_VERTEX;  
       break;  
213      case Tri3Face_Contact:      case Tri3Face_Contact:
       cellType = VTK_LINE;  
       break;  
     case Tri6Face_Contact:  
       cellType = VTK_QUADRATIC_EDGE;  
       break;  
214      case Rec4Face_Contact:      case Rec4Face_Contact:
215        cellType = VTK_LINE;        cellType = VTK_LINE;
       break;  
     case Rec8Face_Contact:  
       cellType = VTK_QUADRATIC_EDGE;  
       break;  
     case Tet4Face_Contact:  
       cellType = VTK_TRIANGLE;  
       break;  
     case Tet10Face_Contact:  
       cellType = VTK_QUADRATIC_TRIANGLE;  
       break;  
     case Hex8Face_Contact:  
       cellType = VTK_QUAD;  
       break;  
     case Hex20Face_Contact:  
       cellType = VTK_QUADRATIC_QUAD;  
       break;  
     default:  
       Finley_ErrorCode=VALUE_ERROR;  
       sprintf(Finley_ErrorMsg,  
           "Element type %s is not supported by VTK",  
           elements->ReferenceElement->Type->Name);  
       return;  
     }  
   
     switch(cellType) {  
     case VTK_VERTEX:  
       numVTKNodesPerElement = 1;  
       strcpy(elemTypeStr, "VTK_VERTEX");  
       break;  
     case VTK_LINE:  
216        numVTKNodesPerElement = 2;        numVTKNodesPerElement = 2;
217        strcpy(elemTypeStr, "VTK_LINE");        strcpy(elemTypeStr, "VTK_LINE");
218        break;        break;
219      case VTK_TRIANGLE:  
220        case Tri3:
221        case Tet4Face:
222        case Tet4Face_Contact:
223          cellType = VTK_TRIANGLE;
224        numVTKNodesPerElement = 3;        numVTKNodesPerElement = 3;
225        strcpy(elemTypeStr, "VTK_TRIANGLE");        strcpy(elemTypeStr, "VTK_TRIANGLE");
226        break;        break;
227      case VTK_QUAD:  
228        case Rec4:
229        case Hex8Face:
230        case Rec4_Contact:
231        case Hex8Face_Contact:
232          cellType = VTK_QUAD;
233        numVTKNodesPerElement = 4;        numVTKNodesPerElement = 4;
234        strcpy(elemTypeStr, "VTK_QUAD");        strcpy(elemTypeStr, "VTK_QUAD");
235        break;        break;
236      case VTK_TETRA:  
237        case Tet4:
238          cellType = VTK_TETRA;
239        numVTKNodesPerElement = 4;        numVTKNodesPerElement = 4;
240        strcpy(elemTypeStr, "VTK_TETRA");        strcpy(elemTypeStr, "VTK_TETRA");
241        break;        break;
242      case VTK_HEXAHEDRON:  
243        case Hex8:
244          cellType = VTK_HEXAHEDRON;
245        numVTKNodesPerElement = 8;        numVTKNodesPerElement = 8;
246        strcpy(elemTypeStr, "VTK_HEXAHEDRON");        strcpy(elemTypeStr, "VTK_HEXAHEDRON");
247        break;        break;
248      case VTK_QUADRATIC_EDGE:  
249        case Line3:
250        case Tri6Face:
251        case Rec8Face:
252        case Line3_Contact:
253        case Tri6Face_Contact:
254        case Rec8Face_Contact:
255          cellType = VTK_QUADRATIC_EDGE;
256        numVTKNodesPerElement = 3;        numVTKNodesPerElement = 3;
257        strcpy(elemTypeStr, "VTK_QUADRATIC_EDGE");        strcpy(elemTypeStr, "VTK_QUADRATIC_EDGE");
258        break;        break;
259      case VTK_QUADRATIC_TRIANGLE:  
260        case Tri6:
261        case Tet10Face:
262        case Tri6_Contact:
263        case Tet10Face_Contact:
264          cellType = VTK_QUADRATIC_TRIANGLE;
265        numVTKNodesPerElement = 6;        numVTKNodesPerElement = 6;
266        strcpy(elemTypeStr, "VTK_QUADRATIC_TRIANGLE");        strcpy(elemTypeStr, "VTK_QUADRATIC_TRIANGLE");
267        break;        break;
268      case VTK_QUADRATIC_QUAD:  
269        case Rec8:
270        case Hex20Face:
271        case Rec8_Contact:
272        case Hex20Face_Contact:
273          cellType = VTK_QUADRATIC_QUAD;
274        numVTKNodesPerElement = 8;        numVTKNodesPerElement = 8;
275        strcpy(elemTypeStr, "VTK_QUADRATIC_QUAD");        strcpy(elemTypeStr, "VTK_QUADRATIC_QUAD");
276        break;        break;
277      case VTK_QUADRATIC_TETRA:  
278        case Tet10:
279          cellType = VTK_QUADRATIC_TETRA;
280        numVTKNodesPerElement = 10;        numVTKNodesPerElement = 10;
281        strcpy(elemTypeStr, "VTK_QUADRATIC_TETRA");        strcpy(elemTypeStr, "VTK_QUADRATIC_TETRA");
282        break;        break;
283      case VTK_QUADRATIC_HEXAHEDRON:  
284        case Hex20:
285          cellType = VTK_QUADRATIC_HEXAHEDRON;
286        numVTKNodesPerElement = 20;        numVTKNodesPerElement = 20;
287        strcpy(elemTypeStr, "VTK_QUADRATIC_HEXAHEDRON");        strcpy(elemTypeStr, "VTK_QUADRATIC_HEXAHEDRON");
288        break;        break;
289      default:  
290        Finley_ErrorCode = VALUE_ERROR;      default:
291        sprintf(Finley_ErrorMsg,        sprintf(error_msg, "saveVTK: Element type %s is not supported by VTK",elements->ReferenceElement->Type->Name);
292            "Cell type %d is not supported by VTK", cellType);        Finley_setError(VALUE_ERROR,error_msg);
293          fclose(fileHandle_p);
294        return;        return;
295      }    }
296      /* xml header */
297      fprintf(fileHandle_p, "<?xml version=\"1.0\"?>\n");
298      fprintf(fileHandle_p, "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\">\n");
299    
300      /* write out the DataArray element for the connectivity */    /* finley uses an unstructured mesh, so UnstructuredGrid *should* work */
301      fprintf(fileHandle_p, "<DataArray "    fprintf(fileHandle_p, "<UnstructuredGrid>\n");
         "Name=\"connectivity\" "  
         "type=\"Int32\" "  
         "format=\"ascii\">\n");  
     int NN = elements->ReferenceElement->Type->numNodes;  
     int counter = 0;  
     for (i = 0; i < numCells; i++) {  
       fprintf(fileHandle_p, "%d ",  
           mesh_p->Elements->Nodes[INDEX2(0, i, NN)]);  
       counter++; /* counter for the number of connectivity points written */  
       /* if the counter gets too big (i.e. the line gets too long),  
        * then add a newline and reset */  
       if (counter > 19) {  
       fprintf(fileHandle_p, "\n");  
       counter = 0;  
       }  
       for (j = 1; j < numVTKNodesPerElement; j++) {  
     fprintf(fileHandle_p,"%d ", mesh_p->Elements->Nodes[INDEX2(j, i, NN)]);  
     counter++;  
     /* if the counter gets too big (i.e. the line gets too long),  
      * then add a newline and reset */  
     if (counter > 19) {  
         fprintf(fileHandle_p, "\n");  
         counter = 0;  
     }  
       }  
     }  
     fprintf(fileHandle_p, "\n");  
     fprintf(fileHandle_p, "</DataArray>\n");  
   
     /* write out the DataArray element for the offsets */  
     fprintf(fileHandle_p, "<DataArray "  
         "Name=\"offsets\" "  
         "type=\"Int32\" "  
         "format=\"ascii\">\n");  
     counter = 0;  /* counter for the number of offsets written to file */  
     for (i=numVTKNodesPerElement; i<=numCells*numVTKNodesPerElement; i+=numVTKNodesPerElement) {  
       fprintf(fileHandle_p, "%d ", i);  
       counter++;  
       /* if the counter gets too big (i.e. the line gets too long),  
        * then add a newline and reset */  
       if (counter > 19) {  
       counter = 0;  
       fprintf(fileHandle_p, "\n");  
       }  
     }  
     fprintf(fileHandle_p, "\n");  
     fprintf(fileHandle_p, "</DataArray>\n");  
   
     /* write out the DataArray element for the types */  
     counter = 0; /* counter for the number of types written to file */  
     fprintf(fileHandle_p, "<DataArray "  
         "Name=\"types\" "  
         "type=\"UInt8\" "  
         "format=\"ascii\">\n");  
     for (i=0; i<numCells; i++) {  
       fprintf(fileHandle_p, "%d ", cellType);  
       counter++;  
       /* if the counter gets too big (i.e. the line gets too long),  
        * then add a newline and reset */  
       if (counter > 30) {  
       counter = 0;  
       fprintf(fileHandle_p, "\n");  
       }  
     }  
     fprintf(fileHandle_p, "\n");  
     fprintf(fileHandle_p, "</DataArray>\n");  
302    
303      /* finish off the <Cells> element */    /* is there only one "piece" to the data?? */
304      fprintf(fileHandle_p, "</Cells>\n");    fprintf(fileHandle_p, "<Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n",numPoints, numCells);
305    }    /* now for the points; equivalent to positions section in saveDX() */
306      /* "The points element explicitly defines coordinates for each point
307       * individually.  It contains one DataArray element describing an array
308       * with three components per value, each specifying the coordinates of one
309       * point" - from Vtk User's Guide
310       */
311      fprintf(fileHandle_p, "<Points>\n");
312      /*
313       * the reason for this if statement is explained in the long comment below
314       */
315      int nDim = mesh_p->Nodes->numDim;
316      fprintf(fileHandle_p, "<DataArray NumberOfComponents=\"%d\" type=\"Float32\" format=\"ascii\">\n",MAX(3,nDim));
317      /* vtk/mayavi doesn't like 2D data, it likes 3D data with a degenerate
318         * third dimension to handle 2D data (like a sheet of paper).  So, if
319         * nDim is 2, we have to append zeros to the array to get this third
320         * dimension, and keep the visualisers happy.
321         * Indeed, if nDim is less than 3, must pad all empty dimensions, so
322         * that the total number of dims is 3.
323      */
324      if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
325         for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
326           if (mesh_p->Nodes->toReduced[i]>=0) {
327              for (j = 0; j < nDim; j++) fprintf(fileHandle_p, " %e",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);
328              for (k=0; k<3-nDim; k++) fprintf(fileHandle_p, " %e",0.);
329              fprintf(fileHandle_p, "\n");
330           }
331         }
332      } else {
333         for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
334           for (j = 0; j < nDim; j++) fprintf(fileHandle_p, " %e",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);
335           for (k=0; k<3-nDim; k++) fprintf(fileHandle_p, " %e",0.);
336           fprintf(fileHandle_p, "\n");
337         }
338      }
339      fprintf(fileHandle_p, "</DataArray>\n");
340      fprintf(fileHandle_p, "</Points>\n");
341    
342    /* data */    /* write out the DataArray element for the connectivity */
343    if (!isEmpty(data_p)) {  
344      int rank = getDataPointRank(data_p);    int NN = elements->ReferenceElement->Type->numNodes;
345      int nComp = getDataPointSize(data_p);    fprintf(fileHandle_p, "<Cells>\n");
346      /* barf if rank is greater than two */    fprintf(fileHandle_p, "<DataArray Name=\"connectivity\" type=\"Int32\" format=\"ascii\">\n");
347      if (rank > 2) {  
348        Finley_ErrorCode = VALUE_ERROR;    if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
349        sprintf(Finley_ErrorMsg,       for (i = 0; i < numCells; i++) {
350            "Vtk can't handle objects with rank greater than 2\n"          for (j = 0; j < numVTKNodesPerElement; j++)
351            "object rank = %d\n", rank);               fprintf(fileHandle_p,"%d ",mesh_p->Nodes->toReduced[elements->Nodes[INDEX2(elements->ReferenceElement->Type->linearNodes[j], i, NN)]]);
352        return;          fprintf(fileHandle_p, "\n");
353      }       }
354      /* if the rank == 0:   --> scalar data    } else if (VTK_QUADRATIC_HEXAHEDRON==cellType) {
355       * if the rank == 1:   --> vector data       for (i = 0; i < numCells; i++) {
356       * if the rank == 2:   --> tensor data            fprintf(fileHandle_p,"%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d\n",
357       */                               elements->Nodes[INDEX2(0, i, NN)],
358      char dataNameStr[31], dataTypeStr[63];                               elements->Nodes[INDEX2(1, i, NN)],
359      int nCompReqd;   /* the number of components required by vtk */                               elements->Nodes[INDEX2(2, i, NN)],
360      if (rank == 0) {                               elements->Nodes[INDEX2(3, i, NN)],
361        strcpy(dataNameStr, "escript_scalar_data");                               elements->Nodes[INDEX2(4, i, NN)],
362        sprintf(dataTypeStr, "Scalars=\"%s\"", dataNameStr);                               elements->Nodes[INDEX2(5, i, NN)],
363        nCompReqd = 1;                               elements->Nodes[INDEX2(6, i, NN)],
364      }                               elements->Nodes[INDEX2(7, i, NN)],
365      else if (rank == 1) {                               elements->Nodes[INDEX2(8, i, NN)],
366        strcpy(dataNameStr, "escript_vector_data");                               elements->Nodes[INDEX2(9, i, NN)],
367        sprintf(dataTypeStr, "Vectors=\"%s\"", dataNameStr);                               elements->Nodes[INDEX2(10, i, NN)],
368        nCompReqd = 3;                               elements->Nodes[INDEX2(11, i, NN)],
369      }                               elements->Nodes[INDEX2(16, i, NN)],
370      else if (rank == 2) {                               elements->Nodes[INDEX2(17, i, NN)],
371        strcpy(dataNameStr, "escript_tensor_data");                               elements->Nodes[INDEX2(18, i, NN)],
372        sprintf(dataTypeStr, "Tensors=\"%s\"", dataNameStr);                               elements->Nodes[INDEX2(19, i, NN)],
373        nCompReqd = 9;                               elements->Nodes[INDEX2(12, i, NN)],
374      }                               elements->Nodes[INDEX2(13, i, NN)],
375      /* if have cell centred data then use CellData element,                               elements->Nodes[INDEX2(14, i, NN)],
376       * if have node centred data, then use PointData element                               elements->Nodes[INDEX2(15, i, NN)]);
377       */       }
378      if (isCellCentered) {    } else if (numVTKNodesPerElement!=NN) {
379        /* now for the cell data */       for (i = 0; i < numCells; i++) {
380        fprintf(fileHandle_p, "<CellData %s>\n", dataTypeStr);          for (j = 0; j < numVTKNodesPerElement; j++) fprintf(fileHandle_p,"%d ", elements->Nodes[INDEX2(elements->ReferenceElement->Type->geoNodes[j], i, NN)]);
381        fprintf(fileHandle_p,          fprintf(fileHandle_p, "\n");
382            "<DataArray "       }
383            "Name=\"%s\" "    } else {
384            "type=\"Float32\" "       for (i = 0; i < numCells; i++) {
385            "NumberOfComponents=\"%d\" "          for (j = 0; j < numVTKNodesPerElement; j++) fprintf(fileHandle_p,"%d ", elements->Nodes[INDEX2(j, i, NN)]);
386            "format=\"ascii\">\n",          fprintf(fileHandle_p, "\n");
387            dataNameStr, nCompReqd);       }
388        int numPointsPerSample = elements->ReferenceElement->numQuadNodes;    }
389        if (numPointsPerSample) {    fprintf(fileHandle_p, "</DataArray>\n");
390      for (i=0; i<numCells; i++) {  
391        values = getSampleData(data_p, i);    /* write out the DataArray element for the offsets */
392        double sampleAvg[nComp];    fprintf(fileHandle_p, "<DataArray Name=\"offsets\" type=\"Int32\" format=\"ascii\">\n");
393        for (k=0; k<nComp; k++) {    for (i=numVTKNodesPerElement; i<=numCells*numVTKNodesPerElement; i+=numVTKNodesPerElement) fprintf(fileHandle_p, "%d\n", i);
394          /* averaging over the number of points in the sample */    fprintf(fileHandle_p, "</DataArray>\n");
         rtmp = 0.;  
         for (j=0; j<numPointsPerSample; j++) {  
           rtmp += values[INDEX2(k,j,nComp)];  
         }  
         sampleAvg[k] = rtmp/numPointsPerSample;  
       }  
       /* if the number of required components is more than the number  
        * of actual components, pad with zeros  
        */  
       /* probably only need to get shape of first element */  
       int shape = getDataPointShape(data_p, 0);  
       if (shape > 3) {  
         Finley_ErrorCode = VALUE_ERROR;  
         sprintf(Finley_ErrorMsg,  
             "shape should be 1, 2, or 3; I got %d\n", shape);  
         return;  
       }  
       /* write the data different ways for scalar, vector and tensor */  
       if (nCompReqd == 1) {  
         fprintf(fileHandle_p, " %f", sampleAvg[0]);  
       }  
       else if (nCompReqd == 3) {  
         int shape = getDataPointShape(data_p, 0);  
         /* write out the data */  
         for (int m=0; m<shape; m++) {  
           fprintf(fileHandle_p, " %f", sampleAvg[m]);  
         }  
         /* pad with zeros */  
         for (int m=0; m<nCompReqd-shape; m++) {  
           fprintf(fileHandle_p, " 0");  
         }  
       }  
       else if (nCompReqd == 9) {  
         /* tensor data, so have a 3x3 matrix to output as a row  
          * of 9 data points */  
         int count = 0;  
         int elems = 0;  
         for (int m=0; m<shape; m++) {  
           for (int n=0; n<shape; n++) {  
         fprintf(fileHandle_p, " %f", sampleAvg[count]);  
         count++;  
         elems++;  
           }  
           for (int n=0; n<3-shape; n++) {  
         fprintf(fileHandle_p, " 0");  
         elems++;  
           }  
         }  
         for (int m=0; m<nCompReqd-elems; m++) {  
           fprintf(fileHandle_p, " 0");  
         }  
       }  
       fprintf(fileHandle_p, "\n");  
     }  
       }  
       fprintf(fileHandle_p, "</DataArray>\n");  
       fprintf(fileHandle_p, "</CellData>\n");  
     } else {  
       /* now for the point data */  
       fprintf(fileHandle_p, "<PointData %s>\n", dataTypeStr);  
       fprintf(fileHandle_p, "<DataArray "  
           "Name=\"%s\" "  
           "type=\"Float32\" "  
           "NumberOfComponents=\"%d\" "  
           "format=\"ascii\">\n",  
           dataNameStr, nCompReqd);  
       for (i=0; i<mesh_p->Nodes->numNodes; i++) {  
     switch (nodetype) {  
     case(FINLEY_DEGREES_OF_FREEDOM):  
       values = getSampleData(data_p,  
                  mesh_p->Nodes->degreeOfFreedom[i]);  
       break;  
     case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):  
       if (mesh_p->Nodes->toReduced[i]>=0) {  
         values = getSampleData(data_p,  
                    mesh_p->Nodes->reducedDegreeOfFreedom[i]);  
       }  
       break;  
     case(FINLEY_NODES):  
       values = getSampleData(data_p,i);  
       break;  
     }  
     /* write out the data */  
     /* if the number of required components is more than the number  
      * of actual components, pad with zeros  
      */  
     /* probably only need to get shape of first element */  
     int shape = getDataPointShape(data_p, 0);  
     if (shape > 3) {  
       Finley_ErrorCode = VALUE_ERROR;  
       sprintf(Finley_ErrorMsg,  
           "shape should be 1, 2, or 3; I got %d\n", shape);  
       return;  
     }  
     /* write the data different ways for scalar, vector and tensor */  
     if (nCompReqd == 1) {  
       fprintf(fileHandle_p, " %f", values[0]);  
     }  
     else if (nCompReqd == 3) {  
       int shape = getDataPointShape(data_p, 0);  
       /* write out the data */  
       for (int m=0; m<shape; m++) {  
         fprintf(fileHandle_p, " %f", values[m]);  
       }  
       /* pad with zeros */  
       for (int m=0; m<nCompReqd-shape; m++) {  
         fprintf(fileHandle_p, " 0");  
       }  
     }  
     else if (nCompReqd == 9) {  
       /* tensor data, so have a 3x3 matrix to output as a row  
        * of 9 data points */  
       int count = 0;  
       int elems = 0;  
       for (int m=0; m<shape; m++) {  
         for (int n=0; n<shape; n++) {  
           fprintf(fileHandle_p, " %f", values[count]);  
           count++;  
           elems++;  
         }  
         for (int n=0; n<3-shape; n++) {  
           fprintf(fileHandle_p, " 0");  
           elems++;  
         }  
       }  
       for (int m=0; m<nCompReqd-elems; m++) {  
         fprintf(fileHandle_p, " 0");  
       }  
     }  
     fprintf(fileHandle_p, "\n");  
       }  
       fprintf(fileHandle_p, "</DataArray>\n");  
       fprintf(fileHandle_p, "</PointData>\n");  
     }  
   }  
395    
396      /* write out the DataArray element for the types */
397      fprintf(fileHandle_p, "<DataArray Name=\"types\" type=\"UInt8\" format=\"ascii\">\n");
398      for (i=0; i<numCells; i++) fprintf(fileHandle_p, "%d\n", cellType);
399      fprintf(fileHandle_p, "</DataArray>\n");
400    
401      /* finish off the <Cells> element */
402      fprintf(fileHandle_p, "</Cells>\n");
403    
404      /* cell data */
405      if (write_celldata) {
406           /* mark the active data arrays */
407           bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
408           fprintf(fileHandle_p, "<CellData");
409           for (i_data =0 ;i_data<num_data;++i_data) {
410                if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data]) {
411                    /* if the rank == 0:   --> scalar data
412                     * if the rank == 1:   --> vector data
413                     * if the rank == 2:   --> tensor data
414                     */
415                    switch(getDataPointRank(data_pp[i_data])) {
416                       case 0:
417                           if (! set_scalar) {
418                                 fprintf(fileHandle_p," Scalars=\"%s\"",names_p[i_data]);
419                                 set_scalar=TRUE;
420                           }
421                           break;
422                       case 1:
423                           if (! set_vector) {
424                                 fprintf(fileHandle_p," Vectors=\"%s\"",names_p[i_data]);
425                                 set_vector=TRUE;
426                           }
427                           break;
428                       case 2:
429                           if (! set_tensor) {
430                                 fprintf(fileHandle_p," Tensors=\"%s\"",names_p[i_data]);
431                                 set_tensor=TRUE;
432                           }
433                           break;
434                       default:
435                           sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);
436                           Finley_setError(VALUE_ERROR,error_msg);
437                           fclose(fileHandle_p);
438                           return;
439                    }
440                }
441           }
442           fprintf(fileHandle_p, ">\n");
443           /* write the arrays */
444           for (i_data =0 ;i_data<num_data;++i_data) {
445              if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data]) {
446                 int numPointsPerSample = elements->ReferenceElement->numQuadNodes;
447                 int rank = getDataPointRank(data_pp[i_data]);
448                 int nComp = getDataPointSize(data_pp[i_data]);
449                 int nCompReqd=1;   /* the number of components required by vtk */
450                 int shape=0;
451                 if (rank == 0) {
452                    nCompReqd = 1;
453                 } else if (rank == 1) {
454                     shape=getDataPointShape(data_pp[i_data], 0);
455                     if  (shape>3) {
456                         Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");
457                         fclose(fileHandle_p);
458                         return;
459                     }
460                     nCompReqd = 3;
461                 } else {
462                     shape=getDataPointShape(data_pp[i_data], 0);
463                     if  (shape>3 || shape != getDataPointShape(data_pp[i_data], 1)) {
464                         Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
465                         fclose(fileHandle_p);
466                         return;
467                     }
468                     nCompReqd = 9;
469                 }
470                 fprintf(fileHandle_p, "<DataArray Name=\"%s\" type=\"Float32\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);
471      
472             double sampleAvg[nComp];
473             for (i=0; i<numCells; i++) {
474                values = getSampleData(data_pp[i_data], i);
475                /* averaging over the number of points in the sample */
476                for (k=0; k<nComp; k++) {
477                   rtmp = 0.;
478                   for (j=0; j<numPointsPerSample; j++) rtmp += values[INDEX2(k,j,nComp)];
479                   sampleAvg[k] = rtmp/numPointsPerSample;
480                 }
481                 /* if the number of required components is more than the number
482                  * of actual components, pad with zeros
483                  */
484                 /* probably only need to get shape of first element */
485                 /* write the data different ways for scalar, vector and tensor */
486                 if (nCompReqd == 1) {
487                   fprintf(fileHandle_p, " %e", sampleAvg[0]);
488                 } else if (nCompReqd == 3) {
489                   /* write out the data */
490                   for (int m=0; m<shape; m++) fprintf(fileHandle_p, " %e", sampleAvg[m]);
491                   for (int m=0; m<nCompReqd-shape; m++) fprintf(fileHandle_p, " %e", 0.);
492                 } else if (nCompReqd == 9) {
493                   /* tensor data, so have a 3x3 matrix to output as a row
494                    * of 9 data points */
495                    int count = 0;
496                    for (int m=0; m<shape; m++) {
497                      for (int n=0; n<shape; n++) {
498                         fprintf(fileHandle_p, " %e", sampleAvg[count]);
499                         count++;
500                      }
501                      for (int n=0; n<3-shape; n++) fprintf(fileHandle_p, " %e", 0.);
502                    }
503                    for (int m=0; m<3-shape; m++)
504                       for (int n=0; n<3; n++) fprintf(fileHandle_p, " %e", 0.);
505                    }
506                  fprintf(fileHandle_p, "\n");
507                 }
508                 fprintf(fileHandle_p, "</DataArray>\n");
509             }
510           }
511           fprintf(fileHandle_p, "</CellData>\n");
512      }
513      /* point data */
514      if (write_pointdata) {
515           /* mark the active data arrays */
516           bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
517           fprintf(fileHandle_p, "<PointData");
518           for (i_data =0 ;i_data<num_data;++i_data) {
519                if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data]) {
520                    /* if the rank == 0:   --> scalar data
521                     * if the rank == 1:   --> vector data
522                     * if the rank == 2:   --> tensor data
523                     */
524                    switch(getDataPointRank(data_pp[i_data])) {
525                       case 0:
526                           if (! set_scalar) {
527                                 fprintf(fileHandle_p," Scalars=\"%s\"",names_p[i_data]);
528                                 set_scalar=TRUE;
529                           }
530                           break;
531                       case 1:
532                           if (! set_vector) {
533                                 fprintf(fileHandle_p," Vectors=\"%s\"",names_p[i_data]);
534                                 set_vector=TRUE;
535                           }
536                           break;
537                       case 2:
538                           if (! set_tensor) {
539                                 fprintf(fileHandle_p," Tensors=\"%s\"",names_p[i_data]);
540                                 set_tensor=TRUE;
541                           }
542                           break;
543                       default:
544                           sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);
545                           Finley_setError(VALUE_ERROR,error_msg);
546                           fclose(fileHandle_p);
547                           return;
548                    }
549                }
550           }
551           fprintf(fileHandle_p, ">\n");
552           /* write the arrays */
553           for (i_data =0 ;i_data<num_data;++i_data) {
554              if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data]) {
555                 int numPointsPerSample = elements->ReferenceElement->numQuadNodes;
556                 int rank = getDataPointRank(data_pp[i_data]);
557                 int nComp = getDataPointSize(data_pp[i_data]);
558                 int nCompReqd=1;   /* the number of components required by vtk */
559                 int shape=0;
560                 if (rank == 0) {
561                    nCompReqd = 1;
562                 } else if (rank == 1) {
563                     shape=getDataPointShape(data_pp[i_data], 0);
564                     if  (shape>3) {
565                         Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");
566                         fclose(fileHandle_p);
567                         return;
568                     }
569                     nCompReqd = 3;
570                 } else {
571                     shape=getDataPointShape(data_pp[i_data], 0);
572                     if  (shape>3 || shape != getDataPointShape(data_pp[i_data], 1)) {
573                         Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
574                         fclose(fileHandle_p);
575                         return;
576                     }
577                     nCompReqd = 9;
578                 }
579                 fprintf(fileHandle_p, "<DataArray Name=\"%s\" type=\"Float32\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);
580             /* write out the data */
581             /* if the number of required components is more than the number
582              * of actual components, pad with zeros
583                   */
584                 bool_t do_write=TRUE;
585                 for (i=0; i<mesh_p->Nodes->numNodes; i++) {
586                   if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
587                    if (mesh_p->Nodes->toReduced[i]>=0) {
588                           switch(getFunctionSpaceType(data_pp[i_data])) {
589                              case FINLEY_DEGREES_OF_FREEDOM:
590                                 values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);
591                                 break;
592                              case FINLEY_REDUCED_DEGREES_OF_FREEDOM:
593                                 values = getSampleData(data_pp[i_data],mesh_p->Nodes->reducedDegreeOfFreedom[i]);
594                                 break;
595                              case FINLEY_NODES:
596                             values = getSampleData(data_pp[i_data],i);
597                                 break;
598                           }
599                           do_write=TRUE;
600                        } else {
601                           do_write=FALSE;
602                        }
603                   } else {
604                        do_write=TRUE;
605                        switch(getFunctionSpaceType(data_pp[i_data])) {
606                           case FINLEY_DEGREES_OF_FREEDOM:
607                          values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]);
608                              break;
609                           case FINLEY_NODES:
610                          values = getSampleData(data_pp[i_data],i);
611                              break;
612                        }
613                   }
614               /* write the data different ways for scalar, vector and tensor */
615                   if (do_write) {
616                  if (nCompReqd == 1) {
617                    fprintf(fileHandle_p, " %e", values[0]);
618                  } else if (nCompReqd == 3) {
619                    for (int m=0; m<shape; m++) fprintf(fileHandle_p, " %e", values[m]);
620                    for (int m=0; m<nCompReqd-shape; m++) fprintf(fileHandle_p, " %e", 0.);
621                  } else if (nCompReqd == 9) {
622                    /* tensor data, so have a 3x3 matrix to output as a row
623                     * of 9 data points */
624                    int count = 0;
625                    for (int m=0; m<shape; m++) {
626                      for (int n=0; n<shape; n++) {
627                        fprintf(fileHandle_p, " %e", values[count]);
628                        count++;
629                      }
630                      for (int n=0; n<3-shape; n++) fprintf(fileHandle_p, " %e", 0.);
631                    }
632                    for (int m=0; m<3-shape; m++)  
633                        for (int n=0; n<3; n++) fprintf(fileHandle_p, " %e", 0.);
634                  }
635                      fprintf(fileHandle_p, "\n");
636                   }
637                 }
638                 fprintf(fileHandle_p, "</DataArray>\n");
639               }
640           }
641           fprintf(fileHandle_p, "</PointData>\n");
642      }
643    /* finish off the piece */    /* finish off the piece */
644    fprintf(fileHandle_p, "</Piece>\n");    fprintf(fileHandle_p, "</Piece>\n");
645    
# Line 614  void Finley_Mesh_saveVTK(const char * fi Line 650  void Finley_Mesh_saveVTK(const char * fi
650    fclose(fileHandle_p);    fclose(fileHandle_p);
651    return;    return;
652  }  }
   
 /*  
  * $Log$  
  * Revision 1.4  2005/05/06 04:26:15  jgs  
  * Merge of development branch back to main trunk on 2005-05-06  
  *  
  * Revision 1.1.2.6  2005/05/06 01:17:19  cochrane  
  * Fixed incorrect reporting of number of components in PointData arrays for  
  * vector data.  
  *  
  * Revision 1.1.2.5  2005/05/05 05:38:44  cochrane  
  * Improved formatting of VTK file output.  
  *  
  * Revision 1.1.2.4  2005/02/22 10:03:54  cochrane  
  * Implementation of writing of vtk xml files from finley.  This function will  
  * require more testing, but on the cases that I have tried (and with the help  
  * of Lutz and mayavi), it looks like it's producing the correct output.  Testing  
  * with more realistic data would be good.  I'm at least confident enough to  
  * commit my changes.  
  *  
  * Revision 1.1.2.3  2005/02/17 05:53:26  gross  
  * some bug in saveDX fixed: in fact the bug was in  
  * DataC/getDataPointShape  
  *  
  * 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.121  
changed lines
  Added in v.153

  ViewVC Help
Powered by ViewVC 1.1.26