/[escript]/branches/domexper/dudley/src/Mesh_saveVTK.c
ViewVC logotype

Diff of /branches/domexper/dudley/src/Mesh_saveVTK.c

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

revision 150 by jgs, Thu Sep 15 03:44:45 2005 UTC revision 153 by jgs, Tue Oct 25 01:51:20 2005 UTC
# Line 29  Line 29 
29    
30  /**************************************************************/  /**************************************************************/
31    
32  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, const dim_t num_data,char* *names_p, escriptDataC* *data_pp) {
33    char error_msg[LenErrorMsg_MAX];    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;
   int i, j, k, numVTKNodesPerElement, isCellCentered=FALSE, nodetype=FINLEY_DEGREES_OF_FREEDOM;  
38    index_t j2;    index_t j2;
39    double* values, rtmp;    double* values, rtmp;
   int nDim = mesh_p->Nodes->numDim;  
   int numPoints=0;  
40        
41    if (nDim==1) {    /* open the file and check handle */
42        Finley_setError(VALUE_ERROR,"saveVTK: 1-dimensional domains are not supported.");  
43        return;    FILE * fileHandle_p = fopen(filename_p, "w");
44      if (fileHandle_p==NULL) {
45        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    /* get a pointer to the relevant mesh component */    /* find the mesh type to be written */
50    if (isEmpty(data_p)) {    int nodetype=FINLEY_DEGREES_OF_FREEDOM;
51      numPoints = mesh_p->Nodes->numNodes;    int elementtype=FINLEY_UNKNOWN;
52      elements=mesh_p->Elements;    bool_t isCellCentered[num_data],write_celldata=FALSE,write_pointdata=FALSE;
53    } else {    for (i_data=0;i_data<num_data;++i_data) {
54      switch(getFunctionSpaceType(data_p)) {       if (! isEmpty(data_pp[i_data])) {
55      case(FINLEY_DEGREES_OF_FREEDOM):          switch(getFunctionSpaceType(data_pp[i_data])) {
56        numPoints = mesh_p->Nodes->numNodes;             case FINLEY_DEGREES_OF_FREEDOM:
57        nodetype = FINLEY_DEGREES_OF_FREEDOM;               nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM;
58        isCellCentered = FALSE;               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
59        elements = mesh_p->Elements;                   elementtype=FINLEY_ELEMENTS;
60        break;               } else {
61      case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):                   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;        numPoints = mesh_p->Nodes->reducedNumNodes;
161        nodetype =FINLEY_REDUCED_DEGREES_OF_FREEDOM;    } else {
       isCellCentered = FALSE;  
       elements = mesh_p->Elements;  
       break;  
     case(FINLEY_NODES):  
       numPoints = mesh_p->Nodes->numNodes;  
       nodetype=FINLEY_NODES;  
       isCellCentered=FALSE;  
       elements=mesh_p->Elements;  
       break;  
     case(FINLEY_ELEMENTS):  
162        numPoints = mesh_p->Nodes->numNodes;        numPoints = mesh_p->Nodes->numNodes;
163        nodetype=FINLEY_NODES;    }
164        isCellCentered=TRUE;    if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;
165      Finley_ElementFile* elements=NULL;
166      switch(elementtype) {
167        case FINLEY_ELEMENTS:
168        elements=mesh_p->Elements;        elements=mesh_p->Elements;
169        break;        break;
170      case(FINLEY_FACE_ELEMENTS):      case FINLEY_FACE_ELEMENTS:
       numPoints = mesh_p->Nodes->numNodes;  
       nodetype=FINLEY_NODES;  
       isCellCentered=TRUE;  
171        elements=mesh_p->FaceElements;        elements=mesh_p->FaceElements;
172        break;        break;
173      case(FINLEY_POINTS):      case FINLEY_POINTS:
       numPoints = mesh_p->Nodes->numNodes;  
       nodetype=FINLEY_NODES;  
       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):  
       numPoints = mesh_p->Nodes->numNodes;  
       nodetype=FINLEY_NODES;  
       isCellCentered=TRUE;  
177        elements=mesh_p->ContactElements;        elements=mesh_p->ContactElements;
178        break;        break;
     default:  
       sprintf(error_msg,"saveVTK: Finley does not know anything about function space type %d",getFunctionSpaceType(data_p));  
       Finley_setError(TYPE_ERROR,error_msg);  
       return;  
     }  
179    }    }
180      if (elements==NULL) {
181    /* the number of points */       Finley_setError(SYSTEM_ERROR,"saveVTK: undefined element file");
182         fclose(fileHandle_p);
183    /* the number of cells */       return;
   if (elements == NULL) {  
     Finley_setError(VALUE_ERROR,"saveVTK: 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    /* open the file and check handle */    int cellType;
188    FILE * fileHandle_p = fopen(filename_p, "w");    ElementTypeId TypeId;
189    if (fileHandle_p==NULL) {    char elemTypeStr[32];
     sprintf(error_msg, "saveVTK: File %s could not be opened for writing.", filename_p);  
     Finley_setError(IO_ERROR,error_msg);  
     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");  
   } else {  
     fprintf(fileHandle_p, "<DataArray "  
         "NumberOfComponents=\"%d\" "  
         "type=\"Float32\" "  
         "format=\"ascii\">\n",  
         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.  
   */  
190    if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {    if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
191       for (i = 0; i < mesh_p->Nodes->numNodes; i++) {       TypeId = elements->LinearReferenceElement->Type->TypeId;
        if (mesh_p->Nodes->toReduced[i]>=0) {  
           for (j = 0; j < nDim; j++)  
             fprintf(fileHandle_p, " %e",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);  
           for (k=0; k<3-nDim; k++) fprintf(fileHandle_p, " %e",0.);  
           fprintf(fileHandle_p, "\n");  
        }  
      }  
192    } else {    } else {
193       for (i = 0; i < mesh_p->Nodes->numNodes; i++) {       TypeId = elements->ReferenceElement->Type->TypeId;
194         for (j = 0; j < nDim; j++)    }
          fprintf(fileHandle_p, " %e",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);  
        for (k=0; k<3-nDim; k++) fprintf(fileHandle_p, " %e",0.);  
        fprintf(fileHandle_p, "\n");  
      }  
   }  
   fprintf(fileHandle_p, "</DataArray>\n");  
   fprintf(fileHandle_p, "</Points>\n");  
   
   /* connections */  
   /* 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 */  
   if (elements!=NULL) {  
     int cellType;  
     ElementTypeId TypeId;  
   
     if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {  
       TypeId = elements->LinearReferenceElement->Type->TypeId;  
     } else {  
       TypeId = elements->ReferenceElement->Type->TypeId;  
     }  
195    
196      switch(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_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_LINE;  
       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:  
       sprintf(error_msg, "saveVTK: Element type %s is not supported by VTK",elements->ReferenceElement->Type->Name);  
       Finley_setError(VALUE_ERROR,error_msg);  
       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        sprintf(error_msg,"saveVTK: Cell type %d is not supported by VTK", cellType);      default:
291          sprintf(error_msg, "saveVTK: Element type %s is not supported by VTK",elements->ReferenceElement->Type->Name);
292        Finley_setError(VALUE_ERROR,error_msg);        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      /* finley uses an unstructured mesh, so UnstructuredGrid *should* work */
301      fprintf(fileHandle_p, "<UnstructuredGrid>\n");
302    
303      /* write out the DataArray element for the connectivity */    /* is there only one "piece" to the data?? */
304      int NN = elements->ReferenceElement->Type->numNodes;    fprintf(fileHandle_p, "<Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n",numPoints, numCells);
305      fprintf(fileHandle_p, "<Cells>\n");    /* now for the points; equivalent to positions section in saveDX() */
306      fprintf(fileHandle_p, "<DataArray "    /* "The points element explicitly defines coordinates for each point
307          "Name=\"connectivity\" "     * individually.  It contains one DataArray element describing an array
308          "type=\"Int32\" "     * with three components per value, each specifying the coordinates of one
309          "format=\"ascii\">\n");     * point" - from Vtk User's Guide
310      if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {     */
311         for (i = 0; i < numCells; i++) {    fprintf(fileHandle_p, "<Points>\n");
312            for (j = 0; j < numVTKNodesPerElement; j++) {    /*
313                 j2=elements->ReferenceElement->Type->linearNodes[j];     * the reason for this if statement is explained in the long comment below
314                 fprintf(fileHandle_p,"%d ",mesh_p->Nodes->toReduced[elements->Nodes[INDEX2(j2, i, NN)]]);     */
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");            fprintf(fileHandle_p, "\n");
330         }         }
     } else if (VTK_QUADRATIC_HEXAHEDRON==cellType) {  
           for (i = 0; i < numCells; i++) {  
             fprintf(fileHandle_p,"%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d\n",  
                                elements->Nodes[INDEX2(0, i, NN)],  
                                elements->Nodes[INDEX2(1, i, NN)],  
                                elements->Nodes[INDEX2(2, i, NN)],  
                                elements->Nodes[INDEX2(3, i, NN)],  
                                elements->Nodes[INDEX2(4, i, NN)],  
                                elements->Nodes[INDEX2(5, i, NN)],  
                                elements->Nodes[INDEX2(6, i, NN)],  
                                elements->Nodes[INDEX2(7, i, NN)],  
                                elements->Nodes[INDEX2(8, i, NN)],  
                                elements->Nodes[INDEX2(9, i, NN)],  
                                elements->Nodes[INDEX2(10, i, NN)],  
                                elements->Nodes[INDEX2(11, i, NN)],  
                                elements->Nodes[INDEX2(16, i, NN)],  
                                elements->Nodes[INDEX2(17, i, NN)],  
                                elements->Nodes[INDEX2(18, i, NN)],  
                                elements->Nodes[INDEX2(19, i, NN)],  
                                elements->Nodes[INDEX2(12, i, NN)],  
                                elements->Nodes[INDEX2(13, i, NN)],  
                                elements->Nodes[INDEX2(14, i, NN)],  
                                elements->Nodes[INDEX2(15, i, NN)]);  
           }  
      } else if (numVTKNodesPerElement!=NN) {  
           for (i = 0; i < numCells; i++) {  
             for (j = 0; j < numVTKNodesPerElement; j++) {  
                  j2=elements->ReferenceElement->Type->geoNodes[j];  
                  fprintf(fileHandle_p,"%d ", elements->Nodes[INDEX2(j2, i, NN)]);  
             }  
             fprintf(fileHandle_p, "\n");  
           }  
      } else {  
           for (i = 0; i < numCells; i++) {  
             for (j = 0; j < numVTKNodesPerElement; j++) fprintf(fileHandle_p,"%d ", elements->Nodes[INDEX2(j, i, NN)]);  
             fprintf(fileHandle_p, "\n");  
           }  
331       }       }
332      fprintf(fileHandle_p, "</DataArray>\n");    } 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      /* write out the DataArray element for the offsets */    /* write out the DataArray element for the connectivity */
     fprintf(fileHandle_p, "<DataArray "  
         "Name=\"offsets\" "  
         "type=\"Int32\" "  
         "format=\"ascii\">\n");  
     for (i=numVTKNodesPerElement; i<=numCells*numVTKNodesPerElement; i+=numVTKNodesPerElement)  
         fprintf(fileHandle_p, "%d\n", i);  
     fprintf(fileHandle_p, "</DataArray>\n");  
   
     /* write out the DataArray element for the types */  
     fprintf(fileHandle_p, "<DataArray "  
         "Name=\"types\" "  
         "type=\"UInt8\" "  
         "format=\"ascii\">\n");  
     for (i=0; i<numCells; i++)  
         fprintf(fileHandle_p, "%d\n", cellType);  
     fprintf(fileHandle_p, "</DataArray>\n");  
343    
344      /* finish off the <Cells> element */    int NN = elements->ReferenceElement->Type->numNodes;
345      fprintf(fileHandle_p, "</Cells>\n");    fprintf(fileHandle_p, "<Cells>\n");
346    }    fprintf(fileHandle_p, "<DataArray Name=\"connectivity\" type=\"Int32\" format=\"ascii\">\n");
347    
348    /* data */    if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
349    if (!isEmpty(data_p)) {       for (i = 0; i < numCells; i++) {
350      int rank = getDataPointRank(data_p);          for (j = 0; j < numVTKNodesPerElement; j++)
351      int nComp = getDataPointSize(data_p);               fprintf(fileHandle_p,"%d ",mesh_p->Nodes->toReduced[elements->Nodes[INDEX2(elements->ReferenceElement->Type->linearNodes[j], i, NN)]]);
352      /* barf if rank is greater than two */          fprintf(fileHandle_p, "\n");
353      if (rank > 2) {       }
354        sprintf(error_msg, "saveVTK: Vtk can't handle objects with rank greater than 2. object rank = %d\n", rank);    } else if (VTK_QUADRATIC_HEXAHEDRON==cellType) {
355        Finley_setError(VALUE_ERROR,error_msg);       for (i = 0; i < numCells; i++) {
356        return;            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      /* if the rank == 0:   --> scalar data                               elements->Nodes[INDEX2(1, i, NN)],
359       * if the rank == 1:   --> vector data                               elements->Nodes[INDEX2(2, i, NN)],
360       * if the rank == 2:   --> tensor data                               elements->Nodes[INDEX2(3, i, NN)],
361       */                               elements->Nodes[INDEX2(4, i, NN)],
362      char dataNameStr[31], dataTypeStr[63];                               elements->Nodes[INDEX2(5, i, NN)],
363      int nCompReqd=1;   /* the number of components required by vtk */                               elements->Nodes[INDEX2(6, i, NN)],
364      if (rank == 0) {                               elements->Nodes[INDEX2(7, i, NN)],
365        strcpy(dataNameStr, "scalar");                               elements->Nodes[INDEX2(8, i, NN)],
366        sprintf(dataTypeStr, "Scalars=\"%s\"", dataNameStr);                               elements->Nodes[INDEX2(9, i, NN)],
367        nCompReqd = 1;                               elements->Nodes[INDEX2(10, i, NN)],
368      }                               elements->Nodes[INDEX2(11, i, NN)],
369      else if (rank == 1) {                               elements->Nodes[INDEX2(16, i, NN)],
370        strcpy(dataNameStr, "vector");                               elements->Nodes[INDEX2(17, i, NN)],
371        sprintf(dataTypeStr, "Vectors=\"%s\"", dataNameStr);                               elements->Nodes[INDEX2(18, i, NN)],
372        nCompReqd = 3;                               elements->Nodes[INDEX2(19, i, NN)],
373      }                               elements->Nodes[INDEX2(12, i, NN)],
374      else if (rank == 2) {                               elements->Nodes[INDEX2(13, i, NN)],
375        strcpy(dataNameStr, "tensor");                               elements->Nodes[INDEX2(14, i, NN)],
376        sprintf(dataTypeStr, "Tensors=\"%s\"", dataNameStr);                               elements->Nodes[INDEX2(15, i, NN)]);
377        nCompReqd = 9;       }
378      }    } else if (numVTKNodesPerElement!=NN) {
379      /* if have cell centred data then use CellData element,       for (i = 0; i < numCells; i++) {
380       * if have node centred data, then use PointData element          for (j = 0; j < numVTKNodesPerElement; j++) fprintf(fileHandle_p,"%d ", elements->Nodes[INDEX2(elements->ReferenceElement->Type->geoNodes[j], i, NN)]);
381       */          fprintf(fileHandle_p, "\n");
382      if (isCellCentered) {       }
383        /* now for the cell data */    } else {
384        fprintf(fileHandle_p, "<CellData %s>\n", dataTypeStr);       for (i = 0; i < numCells; i++) {
385        fprintf(fileHandle_p,          for (j = 0; j < numVTKNodesPerElement; j++) fprintf(fileHandle_p,"%d ", elements->Nodes[INDEX2(j, i, NN)]);
386            "<DataArray "          fprintf(fileHandle_p, "\n");
387            "Name=\"%s\" "       }
388            "type=\"Float32\" "    }
389            "NumberOfComponents=\"%d\" "    fprintf(fileHandle_p, "</DataArray>\n");
390            "format=\"ascii\">\n",  
391            dataNameStr, nCompReqd);    /* write out the DataArray element for the offsets */
392        int numPointsPerSample = elements->ReferenceElement->numQuadNodes;    fprintf(fileHandle_p, "<DataArray Name=\"offsets\" type=\"Int32\" format=\"ascii\">\n");
393        if (numPointsPerSample) {    for (i=numVTKNodesPerElement; i<=numCells*numVTKNodesPerElement; i+=numVTKNodesPerElement) fprintf(fileHandle_p, "%d\n", i);
394      int shape = getDataPointShape(data_p, 0);    fprintf(fileHandle_p, "</DataArray>\n");
395      if (shape > 3) {  
396          sprintf(error_msg,"saveVTK: shape should be 1, 2, or 3; I got %d\n", shape);    /* write out the DataArray element for the types */
397              Finley_setError(VALUE_ERROR,error_msg);    fprintf(fileHandle_p, "<DataArray Name=\"types\" type=\"UInt8\" format=\"ascii\">\n");
398          return;    for (i=0; i<numCells; i++) fprintf(fileHandle_p, "%d\n", cellType);
399      }    fprintf(fileHandle_p, "</DataArray>\n");
400      for (i=0; i<numCells; i++) {  
401        values = getSampleData(data_p, i);    /* finish off the <Cells> element */
402        double sampleAvg[nComp];    fprintf(fileHandle_p, "</Cells>\n");
403        for (k=0; k<nComp; k++) {  
404          /* averaging over the number of points in the sample */    /* cell data */
405          rtmp = 0.;    if (write_celldata) {
406          for (j=0; j<numPointsPerSample; j++)         /* mark the active data arrays */
407            rtmp += values[INDEX2(k,j,nComp)];         bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
408          sampleAvg[k] = rtmp/numPointsPerSample;         fprintf(fileHandle_p, "<CellData");
409        }         for (i_data =0 ;i_data<num_data;++i_data) {
410        /* if the number of required components is more than the number              if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data]) {
411         * of actual components, pad with zeros                  /* if the rank == 0:   --> scalar data
412         */                   * if the rank == 1:   --> vector data
413        /* probably only need to get shape of first element */                   * if the rank == 2:   --> tensor data
414        /* write the data different ways for scalar, vector and tensor */                   */
415        if (nCompReqd == 1) {                  switch(getDataPointRank(data_pp[i_data])) {
416          fprintf(fileHandle_p, " %e", sampleAvg[0]);                     case 0:
417        } else if (nCompReqd == 3) {                         if (! set_scalar) {
418          /* write out the data */                               fprintf(fileHandle_p," Scalars=\"%s\"",names_p[i_data]);
419          for (int m=0; m<shape; m++) {                               set_scalar=TRUE;
420            fprintf(fileHandle_p, " %e", sampleAvg[m]);                         }
421          }                         break;
422          /* pad with zeros */                     case 1:
423          for (int m=0; m<nCompReqd-shape; m++)                         if (! set_vector) {
424            fprintf(fileHandle_p, " %e", 0.);                               fprintf(fileHandle_p," Vectors=\"%s\"",names_p[i_data]);
425        } else if (nCompReqd == 9) {                               set_vector=TRUE;
426          /* tensor data, so have a 3x3 matrix to output as a row                         }
427           * of 9 data points */                         break;
428          int count = 0;                     case 2:
429          for (int m=0; m<shape; m++) {                         if (! set_tensor) {
430            for (int n=0; n<shape; n++) {                               fprintf(fileHandle_p," Tensors=\"%s\"",names_p[i_data]);
431          fprintf(fileHandle_p, " %e", sampleAvg[count]);                               set_tensor=TRUE;
432          count++;                         }
433            }                         break;
434            for (int n=0; n<3-shape; n++)                     default:
435              fprintf(fileHandle_p, " %e", 0.);                         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          for (int m=0; m<3-shape; m++) {                         fclose(fileHandle_p);
438              for (int n=0; n<3; n++)                         return;
439                 fprintf(fileHandle_p, " %e", 0.);                  }
440              }              }
441        }         }
442        fprintf(fileHandle_p, "\n");         fprintf(fileHandle_p, ">\n");
443      }         /* write the arrays */
444        }         for (i_data =0 ;i_data<num_data;++i_data) {
445        fprintf(fileHandle_p, "</DataArray>\n");            if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data]) {
446        fprintf(fileHandle_p, "</CellData>\n");               int numPointsPerSample = elements->ReferenceElement->numQuadNodes;
447      } else {               int rank = getDataPointRank(data_pp[i_data]);
448        /* now for the point data */               int nComp = getDataPointSize(data_pp[i_data]);
449        fprintf(fileHandle_p, "<PointData %s>\n", dataTypeStr);               int nCompReqd=1;   /* the number of components required by vtk */
450        fprintf(fileHandle_p, "<DataArray "               int shape=0;
451            "Name=\"%s\" "               if (rank == 0) {
452            "type=\"Float32\" "                  nCompReqd = 1;
453            "NumberOfComponents=\"%d\" "               } else if (rank == 1) {
454            "format=\"ascii\">\n",                   shape=getDataPointShape(data_pp[i_data], 0);
455            dataNameStr, nCompReqd);                   if  (shape>3) {
456      /* write out the data */                       Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");
457      /* if the number of required components is more than the number                       fclose(fileHandle_p);
458       * of actual components, pad with zeros                       return;
459         */                   }
460        bool_t do_write=TRUE;                   nCompReqd = 3;
461        int shape = getDataPointShape(data_p, 0);               } else {
462        if (shape > 3) {                   shape=getDataPointShape(data_pp[i_data], 0);
463        sprintf(error_msg,"shape should be 1, 2, or 3; I got %d\n", shape);                   if  (shape>3 || shape != getDataPointShape(data_pp[i_data], 1)) {
464            Finley_setError(VALUE_ERROR,error_msg);                       Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
465        return;                       fclose(fileHandle_p);
466        }                       return;
467        for (i=0; i<mesh_p->Nodes->numNodes; i++) {                   }
468      switch (nodetype) {                   nCompReqd = 9;
469      case(FINLEY_DEGREES_OF_FREEDOM):               }
470        values = getSampleData(data_p,mesh_p->Nodes->degreeOfFreedom[i]);               fprintf(fileHandle_p, "<DataArray Name=\"%s\" type=\"Float32\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd);
471        break;    
472      case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):           double sampleAvg[nComp];
473        if (mesh_p->Nodes->toReduced[i]>=0) {           for (i=0; i<numCells; i++) {
474              do_write=TRUE;              values = getSampleData(data_pp[i_data], i);
475          values = getSampleData(data_p,mesh_p->Nodes->reducedDegreeOfFreedom[i]);              /* averaging over the number of points in the sample */
476        } else {              for (k=0; k<nComp; k++) {
477              do_write=FALSE;                 rtmp = 0.;
478            }                 for (j=0; j<numPointsPerSample; j++) rtmp += values[INDEX2(k,j,nComp)];
479        break;                 sampleAvg[k] = rtmp/numPointsPerSample;
480      case(FINLEY_NODES):               }
481        values = getSampleData(data_p,i);               /* if the number of required components is more than the number
482        break;                * of actual components, pad with zeros
483      }                */
484      /* write the data different ways for scalar, vector and tensor */               /* probably only need to get shape of first element */
485          if (do_write) {               /* write the data different ways for scalar, vector and tensor */
486          if (nCompReqd == 1) {               if (nCompReqd == 1) {
487            fprintf(fileHandle_p, " %e", values[0]);                 fprintf(fileHandle_p, " %e", sampleAvg[0]);
488          }               } else if (nCompReqd == 3) {
489          else if (nCompReqd == 3) {                 /* write out the data */
490            /* write out the data */                 for (int m=0; m<shape; m++) fprintf(fileHandle_p, " %e", sampleAvg[m]);
491            for (int m=0; m<shape; m++)                 for (int m=0; m<nCompReqd-shape; m++) fprintf(fileHandle_p, " %e", 0.);
492              fprintf(fileHandle_p, " %e", values[m]);               } else if (nCompReqd == 9) {
493            /* pad with zeros */                 /* tensor data, so have a 3x3 matrix to output as a row
494            for (int m=0; m<nCompReqd-shape; m++)                  * of 9 data points */
495                fprintf(fileHandle_p, " %e", 0.);                  int count = 0;
496          }                  for (int m=0; m<shape; m++) {
497          else if (nCompReqd == 9) {                    for (int n=0; n<shape; n++) {
498            /* tensor data, so have a 3x3 matrix to output as a row                       fprintf(fileHandle_p, " %e", sampleAvg[count]);
499             * of 9 data points */                       count++;
500            int count = 0;                    }
501            for (int m=0; m<shape; m++) {                    for (int n=0; n<3-shape; n++) fprintf(fileHandle_p, " %e", 0.);
502              for (int n=0; n<shape; n++) {                  }
503                fprintf(fileHandle_p, " %e", values[count]);                  for (int m=0; m<3-shape; m++)
504                count++;                     for (int n=0; n<3; n++) fprintf(fileHandle_p, " %e", 0.);
505              }                  }
506              for (int n=0; n<3-shape; n++)                fprintf(fileHandle_p, "\n");
507                fprintf(fileHandle_p, " %e", 0.);               }
508            }               fprintf(fileHandle_p, "</DataArray>\n");
           for (int m=0; m<3-shape; m++)  {  
               for (int n=0; n<3; n++)  
                   fprintf(fileHandle_p, " %e", 0.);  
               }  
         }  
         fprintf(fileHandle_p, "\n");  
509           }           }
510        }         }
511        fprintf(fileHandle_p, "</DataArray>\n");         fprintf(fileHandle_p, "</CellData>\n");
512        fprintf(fileHandle_p, "</PointData>\n");    }
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 630  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  }  }
   
 /*  
  * Revision 1.6  2005/08/12 01:45:43  jgs  
  * erge of development branch dev-02 back to main trunk on 2005-08-12  
  *  
  * Revision 1.5.2.4  2005/09/09 08:15:17  gross  
  * bugs in vtk and dx writer fixed  
  *  
  * Revision 1.5.2.3  2005/09/08 08:28:39  gross  
  * some cleanup in savevtk  
  *  
  * Revision 1.5.2.2  2005/09/07 06:26:20  gross  
  * the solver from finley are put into the standalone package paso now  
  *  
  * Revision 1.5.2.1  2005/08/10 06:14:37  gross  
  * QUADRATIC HEXAHEDRON elements fixed  
  *  
  * Revision 1.5  2005/07/08 04:07:55  jgs  
  * Merge of development branch back to main trunk on 2005-07-08  
  *  
  * 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.7  2005/06/29 02:34:54  gross  
  * some changes towards 64 integers in finley  
  *  
  * 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.150  
changed lines
  Added in v.153

  ViewVC Help
Powered by ViewVC 1.1.26