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

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

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

revision 147 by jgs, Fri Aug 12 01:45:47 2005 UTC revision 150 by jgs, Thu Sep 15 03:44:45 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    /**************************************************************/
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, escriptDataC* data_p) {
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;    Finley_ElementFile* elements=NULL;
37    char elemTypeStr[32];    char elemTypeStr[32];
38    int i, j, k, numVTKNodesPerElement, isCellCentered=FALSE, nodetype=FINLEY_DEGREES_OF_FREEDOM;    int i, j, k, numVTKNodesPerElement, isCellCentered=FALSE, nodetype=FINLEY_DEGREES_OF_FREEDOM;
39      index_t j2;
40    double* values, rtmp;    double* values, rtmp;
41    int nDim = mesh_p->Nodes->numDim;    int nDim = mesh_p->Nodes->numDim;
42      int numPoints=0;
43      
44      if (nDim==1) {
45          Finley_setError(VALUE_ERROR,"saveVTK: 1-dimensional domains are not supported.");
46          return;
47      }
48    /* get a pointer to the relevant mesh component */    /* get a pointer to the relevant mesh component */
49    if (isEmpty(data_p)) {    if (isEmpty(data_p)) {
50        numPoints = mesh_p->Nodes->numNodes;
51      elements=mesh_p->Elements;      elements=mesh_p->Elements;
52    } else {    } else {
53      switch(getFunctionSpaceType(data_p)) {      switch(getFunctionSpaceType(data_p)) {
54      case(FINLEY_DEGREES_OF_FREEDOM):      case(FINLEY_DEGREES_OF_FREEDOM):
55          numPoints = mesh_p->Nodes->numNodes;
56        nodetype = FINLEY_DEGREES_OF_FREEDOM;        nodetype = FINLEY_DEGREES_OF_FREEDOM;
57        isCellCentered = FALSE;        isCellCentered = FALSE;
58        elements = mesh_p->Elements;        elements = mesh_p->Elements;
59        break;        break;
60      case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):      case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):
61        Finley_ErrorCode=VALUE_ERROR;        numPoints = mesh_p->Nodes->reducedNumNodes;
62        sprintf(Finley_ErrorMsg,        nodetype =FINLEY_REDUCED_DEGREES_OF_FREEDOM;
63            "Reduced degrees of freedom is not yet "        isCellCentered = FALSE;
64            "implemented for saving vtk files\n");        elements = mesh_p->Elements;
65        return;        break;
66      case(FINLEY_NODES):      case(FINLEY_NODES):
67          numPoints = mesh_p->Nodes->numNodes;
68        nodetype=FINLEY_NODES;        nodetype=FINLEY_NODES;
69        isCellCentered=FALSE;        isCellCentered=FALSE;
70        elements=mesh_p->Elements;        elements=mesh_p->Elements;
71        break;        break;
72      case(FINLEY_ELEMENTS):      case(FINLEY_ELEMENTS):
73          numPoints = mesh_p->Nodes->numNodes;
74          nodetype=FINLEY_NODES;
75        isCellCentered=TRUE;        isCellCentered=TRUE;
76        elements=mesh_p->Elements;        elements=mesh_p->Elements;
77        break;        break;
78      case(FINLEY_FACE_ELEMENTS):      case(FINLEY_FACE_ELEMENTS):
79          numPoints = mesh_p->Nodes->numNodes;
80          nodetype=FINLEY_NODES;
81        isCellCentered=TRUE;        isCellCentered=TRUE;
82        elements=mesh_p->FaceElements;        elements=mesh_p->FaceElements;
83        break;        break;
84      case(FINLEY_POINTS):      case(FINLEY_POINTS):
85          numPoints = mesh_p->Nodes->numNodes;
86          nodetype=FINLEY_NODES;
87        isCellCentered=TRUE;        isCellCentered=TRUE;
88        elements=mesh_p->Points;        elements=mesh_p->Points;
89        break;        break;
90      case(FINLEY_CONTACT_ELEMENTS_1):      case(FINLEY_CONTACT_ELEMENTS_1):
91      case(FINLEY_CONTACT_ELEMENTS_2):      case(FINLEY_CONTACT_ELEMENTS_2):
92          numPoints = mesh_p->Nodes->numNodes;
93          nodetype=FINLEY_NODES;
94        isCellCentered=TRUE;        isCellCentered=TRUE;
95        elements=mesh_p->ContactElements;        elements=mesh_p->ContactElements;
96        break;        break;
97      default:      default:
98        Finley_ErrorCode=TYPE_ERROR;        sprintf(error_msg,"saveVTK: Finley does not know anything about function space type %d",getFunctionSpaceType(data_p));
99        sprintf(Finley_ErrorMsg,        Finley_setError(TYPE_ERROR,error_msg);
           "Finley does not know anything about function space type %d",  
           getFunctionSpaceType(data_p));  
100        return;        return;
101      }      }
102    }    }
103    
104    /* the number of points */    /* the number of points */
   int numPoints = mesh_p->Nodes->numNodes;  
105    
106    /* the number of cells */    /* the number of cells */
107    if (elements == NULL) {    if (elements == NULL) {
108      Finley_ErrorCode = VALUE_ERROR;      Finley_setError(VALUE_ERROR,"saveVTK: elements object is NULL; cannot proceed");
     sprintf(Finley_ErrorMsg,  
         "elements object is NULL; cannot proceed");  
109      return;      return;
110    }    }
111    int numCells = elements->numElements;      int numCells = elements->numElements;  
     
112    /* open the file and check handle */    /* open the file and check handle */
113    FILE * fileHandle_p = fopen(filename_p, "w");    FILE * fileHandle_p = fopen(filename_p, "w");
114    if (fileHandle_p==NULL) {    if (fileHandle_p==NULL) {
115      Finley_ErrorCode=IO_ERROR;      sprintf(error_msg, "saveVTK: File %s could not be opened for writing.", filename_p);
116      sprintf(Finley_ErrorMsg,      Finley_setError(IO_ERROR,error_msg);
         "File %s could not be opened for writing.", filename_p);  
117      return;      return;
118    }    }
119    /* xml header */    /* xml header */
120    fprintf(fileHandle_p, "<?xml version=\"1.0\"?>\n");    fprintf(fileHandle_p, "<?xml version=\"1.0\"?>\n");
121    fprintf(fileHandle_p,    fprintf(fileHandle_p, "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\">\n");
       "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\">\n");  
122    
123    /* finley uses an unstructured mesh, so UnstructuredGrid *should* work */    /* finley uses an unstructured mesh, so UnstructuredGrid *should* work */
124    fprintf(fileHandle_p, "<UnstructuredGrid>\n");    fprintf(fileHandle_p, "<UnstructuredGrid>\n");
125    
126    /* is there only one "piece" to the data?? */    /* is there only one "piece" to the data?? */
127    fprintf(fileHandle_p, "<Piece "    fprintf(fileHandle_p, "<Piece "
128        "NumberOfPoints=\"%d\" "            "NumberOfPoints=\"%d\" "
129        "NumberOfCells=\"%d\">\n",            "NumberOfCells=\"%d\">\n",
130        numPoints, numCells);            numPoints, numCells);
   
131    /* now for the points; equivalent to positions section in saveDX() */    /* now for the points; equivalent to positions section in saveDX() */
132    /* "The points element explicitly defines coordinates for each point    /* "The points element explicitly defines coordinates for each point
133     * individually.  It contains one DataArray element describing an array     * individually.  It contains one DataArray element describing an array
# Line 129  void Finley_Mesh_saveVTK(const char * fi Line 150  void Finley_Mesh_saveVTK(const char * fi
150          "format=\"ascii\">\n",          "format=\"ascii\">\n",
151          nDim);          nDim);
152    }    }
153    for (i = 0; i < mesh_p->Nodes->numNodes; i++) {    /* vtk/mayavi doesn't like 2D data, it likes 3D data with a degenerate
154      fprintf(fileHandle_p,       * third dimension to handle 2D data (like a sheet of paper).  So, if
155          "%e", mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)]);       * nDim is 2, we have to append zeros to the array to get this third
156      for (j = 1; j < nDim; j++) {       * dimension, and keep the visualisers happy.
157        fprintf(fileHandle_p,       * Indeed, if nDim is less than 3, must pad all empty dimensions, so
158            " %f",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);       * that the total number of dims is 3.
159        /* vtk/mayavi doesn't like 2D data, it likes 3D data with a degenerate    */
160         * third dimension to handle 2D data (like a sheet of paper).  So, if    if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
161         * nDim is 2, we have to append zeros to the array to get this third       for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
162         * dimension, and keep the visualisers happy.         if (mesh_p->Nodes->toReduced[i]>=0) {
163         * Indeed, if nDim is less than 3, must pad all empty dimensions, so            for (j = 0; j < nDim; j++)
164         * that the total number of dims is 3.              fprintf(fileHandle_p, " %e",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);
165         */            for (k=0; k<3-nDim; k++) fprintf(fileHandle_p, " %e",0.);
166        if (nDim < 3) {            fprintf(fileHandle_p, "\n");
167      for (k=0; k<3-nDim; k++) {         }
168        fprintf(fileHandle_p, " 0");       }
169      }    } else {
170        }       for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
171      }         for (j = 0; j < nDim; j++)
172      fprintf(fileHandle_p, "\n");           fprintf(fileHandle_p, " %e",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);
173           for (k=0; k<3-nDim; k++) fprintf(fileHandle_p, " %e",0.);
174           fprintf(fileHandle_p, "\n");
175         }
176    }    }
177    fprintf(fileHandle_p, "</DataArray>\n");    fprintf(fileHandle_p, "</DataArray>\n");
178    fprintf(fileHandle_p, "</Points>\n");    fprintf(fileHandle_p, "</Points>\n");
# Line 163  void Finley_Mesh_saveVTK(const char * fi Line 187  void Finley_Mesh_saveVTK(const char * fi
187     * specifies the type of each cell.     * specifies the type of each cell.
188     */     */
189    /* if no element table is present jump over the connection table */    /* if no element table is present jump over the connection table */
   int cellType;  
190    if (elements!=NULL) {    if (elements!=NULL) {
191      fprintf(fileHandle_p, "<Cells>\n");      int cellType;
192      ElementTypeId TypeId = elements->ReferenceElement->Type->TypeId;      ElementTypeId TypeId;
193    
194        if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
195          TypeId = elements->LinearReferenceElement->Type->TypeId;
196        } else {
197          TypeId = elements->ReferenceElement->Type->TypeId;
198        }
199    
200      switch(TypeId) {      switch(TypeId) {
201      case Point1:      case Point1:
202        cellType = VTK_VERTEX;        cellType = VTK_VERTEX;
# Line 226  void Finley_Mesh_saveVTK(const char * fi Line 256  void Finley_Mesh_saveVTK(const char * fi
256        cellType = VTK_QUADRATIC_TRIANGLE;        cellType = VTK_QUADRATIC_TRIANGLE;
257        break;        break;
258      case Hex8Face:      case Hex8Face:
259        cellType = VTK_QUADRATIC_QUAD;        cellType = VTK_QUAD;
260        break;        break;
261      case Hex20Face:      case Hex20Face:
262        cellType = VTK_QUADRATIC_QUAD;        cellType = VTK_QUADRATIC_QUAD;
# Line 241  void Finley_Mesh_saveVTK(const char * fi Line 271  void Finley_Mesh_saveVTK(const char * fi
271        cellType = VTK_QUADRATIC_EDGE;        cellType = VTK_QUADRATIC_EDGE;
272        break;        break;
273      case Tri3_Contact:      case Tri3_Contact:
274        cellType = VTK_TRIANGLE;        cellType = VTK_LINE;
275        break;        break;
276      case Tri6_Contact:      case Tri6_Contact:
277        cellType = VTK_QUADRATIC_TRIANGLE;        cellType = VTK_QUADRATIC_TRIANGLE;
# Line 283  void Finley_Mesh_saveVTK(const char * fi Line 313  void Finley_Mesh_saveVTK(const char * fi
313        cellType = VTK_QUADRATIC_QUAD;        cellType = VTK_QUADRATIC_QUAD;
314        break;        break;
315      default:      default:
316        Finley_ErrorCode=VALUE_ERROR;        sprintf(error_msg, "saveVTK: Element type %s is not supported by VTK",elements->ReferenceElement->Type->Name);
317        sprintf(Finley_ErrorMsg,        Finley_setError(VALUE_ERROR,error_msg);
           "Element type %s is not supported by VTK",  
           elements->ReferenceElement->Type->Name);  
318        return;        return;
319      }      }
320    
# Line 336  void Finley_Mesh_saveVTK(const char * fi Line 364  void Finley_Mesh_saveVTK(const char * fi
364        strcpy(elemTypeStr, "VTK_QUADRATIC_HEXAHEDRON");        strcpy(elemTypeStr, "VTK_QUADRATIC_HEXAHEDRON");
365        break;        break;
366      default:      default:
367        Finley_ErrorCode = VALUE_ERROR;        sprintf(error_msg,"saveVTK: Cell type %d is not supported by VTK", cellType);
368        sprintf(Finley_ErrorMsg,        Finley_setError(VALUE_ERROR,error_msg);
           "Cell type %d is not supported by VTK", cellType);  
369        return;        return;
370      }      }
371    
372      /* write out the DataArray element for the connectivity */      /* write out the DataArray element for the connectivity */
373        int NN = elements->ReferenceElement->Type->numNodes;
374        fprintf(fileHandle_p, "<Cells>\n");
375      fprintf(fileHandle_p, "<DataArray "      fprintf(fileHandle_p, "<DataArray "
376          "Name=\"connectivity\" "          "Name=\"connectivity\" "
377          "type=\"Int32\" "          "type=\"Int32\" "
378          "format=\"ascii\">\n");          "format=\"ascii\">\n");
379      int NN = elements->ReferenceElement->Type->numNodes;      if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
     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",  
                             mesh_p->Elements->Nodes[INDEX2(0, i, NN)],  
                             mesh_p->Elements->Nodes[INDEX2(1, i, NN)],  
                             mesh_p->Elements->Nodes[INDEX2(2, i, NN)],  
                             mesh_p->Elements->Nodes[INDEX2(3, i, NN)],  
                             mesh_p->Elements->Nodes[INDEX2(4, i, NN)],  
                             mesh_p->Elements->Nodes[INDEX2(5, i, NN)],  
                             mesh_p->Elements->Nodes[INDEX2(6, i, NN)],  
                             mesh_p->Elements->Nodes[INDEX2(7, i, NN)],  
                             mesh_p->Elements->Nodes[INDEX2(8, i, NN)],  
                             mesh_p->Elements->Nodes[INDEX2(9, i, NN)],  
                             mesh_p->Elements->Nodes[INDEX2(10, i, NN)],  
                             mesh_p->Elements->Nodes[INDEX2(11, i, NN)],  
                             mesh_p->Elements->Nodes[INDEX2(16, i, NN)],  
                             mesh_p->Elements->Nodes[INDEX2(17, i, NN)],  
                             mesh_p->Elements->Nodes[INDEX2(18, i, NN)],  
                             mesh_p->Elements->Nodes[INDEX2(19, i, NN)],  
                             mesh_p->Elements->Nodes[INDEX2(12, i, NN)],  
                             mesh_p->Elements->Nodes[INDEX2(13, i, NN)],  
                             mesh_p->Elements->Nodes[INDEX2(14, i, NN)],  
                             mesh_p->Elements->Nodes[INDEX2(15, i, NN)]);  
          fprintf(fileHandle_p, "\n");  
        }  
     } else {  
380         for (i = 0; i < numCells; i++) {         for (i = 0; i < numCells; i++) {
381           for (j = 0; j < numVTKNodesPerElement; j++) fprintf(fileHandle_p,"%d ", mesh_p->Elements->Nodes[INDEX2(j, i, NN)]);            for (j = 0; j < numVTKNodesPerElement; j++) {
382           fprintf(fileHandle_p, "\n");                 j2=elements->ReferenceElement->Type->linearNodes[j];
383                   fprintf(fileHandle_p,"%d ",mesh_p->Nodes->toReduced[elements->Nodes[INDEX2(j2, i, NN)]]);
384              }
385              fprintf(fileHandle_p, "\n");
386         }         }
387      }      } else if (VTK_QUADRATIC_HEXAHEDRON==cellType) {
388      fprintf(fileHandle_p, "\n");            for (i = 0; i < numCells; i++) {
389                fprintf(fileHandle_p,"%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d\n",
390                                   elements->Nodes[INDEX2(0, i, NN)],
391                                   elements->Nodes[INDEX2(1, i, NN)],
392                                   elements->Nodes[INDEX2(2, i, NN)],
393                                   elements->Nodes[INDEX2(3, i, NN)],
394                                   elements->Nodes[INDEX2(4, i, NN)],
395                                   elements->Nodes[INDEX2(5, i, NN)],
396                                   elements->Nodes[INDEX2(6, i, NN)],
397                                   elements->Nodes[INDEX2(7, i, NN)],
398                                   elements->Nodes[INDEX2(8, i, NN)],
399                                   elements->Nodes[INDEX2(9, i, NN)],
400                                   elements->Nodes[INDEX2(10, i, NN)],
401                                   elements->Nodes[INDEX2(11, i, NN)],
402                                   elements->Nodes[INDEX2(16, i, NN)],
403                                   elements->Nodes[INDEX2(17, i, NN)],
404                                   elements->Nodes[INDEX2(18, i, NN)],
405                                   elements->Nodes[INDEX2(19, i, NN)],
406                                   elements->Nodes[INDEX2(12, i, NN)],
407                                   elements->Nodes[INDEX2(13, i, NN)],
408                                   elements->Nodes[INDEX2(14, i, NN)],
409                                   elements->Nodes[INDEX2(15, i, NN)]);
410              }
411         } else if (numVTKNodesPerElement!=NN) {
412              for (i = 0; i < numCells; i++) {
413                for (j = 0; j < numVTKNodesPerElement; j++) {
414                     j2=elements->ReferenceElement->Type->geoNodes[j];
415                     fprintf(fileHandle_p,"%d ", elements->Nodes[INDEX2(j2, i, NN)]);
416                }
417                fprintf(fileHandle_p, "\n");
418              }
419         } else {
420              for (i = 0; i < numCells; i++) {
421                for (j = 0; j < numVTKNodesPerElement; j++) fprintf(fileHandle_p,"%d ", elements->Nodes[INDEX2(j, i, NN)]);
422                fprintf(fileHandle_p, "\n");
423              }
424         }
425      fprintf(fileHandle_p, "</DataArray>\n");      fprintf(fileHandle_p, "</DataArray>\n");
426    
427      /* write out the DataArray element for the offsets */      /* write out the DataArray element for the offsets */
# Line 387  void Finley_Mesh_saveVTK(const char * fi Line 429  void Finley_Mesh_saveVTK(const char * fi
429          "Name=\"offsets\" "          "Name=\"offsets\" "
430          "type=\"Int32\" "          "type=\"Int32\" "
431          "format=\"ascii\">\n");          "format=\"ascii\">\n");
432      int counter = 0;  /* counter for the number of offsets written to file */      for (i=numVTKNodesPerElement; i<=numCells*numVTKNodesPerElement; i+=numVTKNodesPerElement)
433      for (i=numVTKNodesPerElement; i<=numCells*numVTKNodesPerElement; i+=numVTKNodesPerElement) {          fprintf(fileHandle_p, "%d\n", i);
       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");  
434      fprintf(fileHandle_p, "</DataArray>\n");      fprintf(fileHandle_p, "</DataArray>\n");
435    
436      /* write out the DataArray element for the types */      /* write out the DataArray element for the types */
     counter = 0; /* counter for the number of types written to file */  
437      fprintf(fileHandle_p, "<DataArray "      fprintf(fileHandle_p, "<DataArray "
438          "Name=\"types\" "          "Name=\"types\" "
439          "type=\"UInt8\" "          "type=\"UInt8\" "
440          "format=\"ascii\">\n");          "format=\"ascii\">\n");
441      for (i=0; i<numCells; i++) {      for (i=0; i<numCells; i++)
442        fprintf(fileHandle_p, "%d ", cellType);          fprintf(fileHandle_p, "%d\n", 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");  
443      fprintf(fileHandle_p, "</DataArray>\n");      fprintf(fileHandle_p, "</DataArray>\n");
444    
445      /* finish off the <Cells> element */      /* finish off the <Cells> element */
# Line 430  void Finley_Mesh_saveVTK(const char * fi Line 452  void Finley_Mesh_saveVTK(const char * fi
452      int nComp = getDataPointSize(data_p);      int nComp = getDataPointSize(data_p);
453      /* barf if rank is greater than two */      /* barf if rank is greater than two */
454      if (rank > 2) {      if (rank > 2) {
455        Finley_ErrorCode = VALUE_ERROR;        sprintf(error_msg, "saveVTK: Vtk can't handle objects with rank greater than 2. object rank = %d\n", rank);
456        sprintf(Finley_ErrorMsg,        Finley_setError(VALUE_ERROR,error_msg);
           "Vtk can't handle objects with rank greater than 2\n"  
           "object rank = %d\n", rank);  
457        return;        return;
458      }      }
459      /* if the rank == 0:   --> scalar data      /* if the rank == 0:   --> scalar data
# Line 443  void Finley_Mesh_saveVTK(const char * fi Line 463  void Finley_Mesh_saveVTK(const char * fi
463      char dataNameStr[31], dataTypeStr[63];      char dataNameStr[31], dataTypeStr[63];
464      int nCompReqd=1;   /* the number of components required by vtk */      int nCompReqd=1;   /* the number of components required by vtk */
465      if (rank == 0) {      if (rank == 0) {
466        strcpy(dataNameStr, "escript_scalar_data");        strcpy(dataNameStr, "scalar");
467        sprintf(dataTypeStr, "Scalars=\"%s\"", dataNameStr);        sprintf(dataTypeStr, "Scalars=\"%s\"", dataNameStr);
468        nCompReqd = 1;        nCompReqd = 1;
469      }      }
470      else if (rank == 1) {      else if (rank == 1) {
471        strcpy(dataNameStr, "escript_vector_data");        strcpy(dataNameStr, "vector");
472        sprintf(dataTypeStr, "Vectors=\"%s\"", dataNameStr);        sprintf(dataTypeStr, "Vectors=\"%s\"", dataNameStr);
473        nCompReqd = 3;        nCompReqd = 3;
474      }      }
475      else if (rank == 2) {      else if (rank == 2) {
476        strcpy(dataNameStr, "escript_tensor_data");        strcpy(dataNameStr, "tensor");
477        sprintf(dataTypeStr, "Tensors=\"%s\"", dataNameStr);        sprintf(dataTypeStr, "Tensors=\"%s\"", dataNameStr);
478        nCompReqd = 9;        nCompReqd = 9;
479      }      }
# Line 472  void Finley_Mesh_saveVTK(const char * fi Line 492  void Finley_Mesh_saveVTK(const char * fi
492            dataNameStr, nCompReqd);            dataNameStr, nCompReqd);
493        int numPointsPerSample = elements->ReferenceElement->numQuadNodes;        int numPointsPerSample = elements->ReferenceElement->numQuadNodes;
494        if (numPointsPerSample) {        if (numPointsPerSample) {
495        int shape = getDataPointShape(data_p, 0);
496        if (shape > 3) {
497            sprintf(error_msg,"saveVTK: shape should be 1, 2, or 3; I got %d\n", shape);
498                Finley_setError(VALUE_ERROR,error_msg);
499            return;
500        }
501      for (i=0; i<numCells; i++) {      for (i=0; i<numCells; i++) {
502        values = getSampleData(data_p, i);        values = getSampleData(data_p, i);
503        double sampleAvg[nComp];        double sampleAvg[nComp];
504        for (k=0; k<nComp; k++) {        for (k=0; k<nComp; k++) {
505          /* averaging over the number of points in the sample */          /* averaging over the number of points in the sample */
506          rtmp = 0.;          rtmp = 0.;
507          for (j=0; j<numPointsPerSample; j++) {          for (j=0; j<numPointsPerSample; j++)
508            rtmp += values[INDEX2(k,j,nComp)];            rtmp += values[INDEX2(k,j,nComp)];
         }  
509          sampleAvg[k] = rtmp/numPointsPerSample;          sampleAvg[k] = rtmp/numPointsPerSample;
510        }        }
511        /* if the number of required components is more than the number        /* if the number of required components is more than the number
512         * of actual components, pad with zeros         * of actual components, pad with zeros
513         */         */
514        /* probably only need to get shape of first element */        /* 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;  
       }  
515        /* write the data different ways for scalar, vector and tensor */        /* write the data different ways for scalar, vector and tensor */
516        if (nCompReqd == 1) {        if (nCompReqd == 1) {
517          fprintf(fileHandle_p, " %f", sampleAvg[0]);          fprintf(fileHandle_p, " %e", sampleAvg[0]);
518        }        } else if (nCompReqd == 3) {
       else if (nCompReqd == 3) {  
         int shape = getDataPointShape(data_p, 0);  
519          /* write out the data */          /* write out the data */
520          for (int m=0; m<shape; m++) {          for (int m=0; m<shape; m++) {
521            fprintf(fileHandle_p, " %f", sampleAvg[m]);            fprintf(fileHandle_p, " %e", sampleAvg[m]);
522          }          }
523          /* pad with zeros */          /* pad with zeros */
524          for (int m=0; m<nCompReqd-shape; m++) {          for (int m=0; m<nCompReqd-shape; m++)
525            fprintf(fileHandle_p, " 0");            fprintf(fileHandle_p, " %e", 0.);
526          }        } else if (nCompReqd == 9) {
       }  
       else if (nCompReqd == 9) {  
527          /* tensor data, so have a 3x3 matrix to output as a row          /* tensor data, so have a 3x3 matrix to output as a row
528           * of 9 data points */           * of 9 data points */
529          int count = 0;          int count = 0;
         int elems = 0;  
530          for (int m=0; m<shape; m++) {          for (int m=0; m<shape; m++) {
531            for (int n=0; n<shape; n++) {            for (int n=0; n<shape; n++) {
532          fprintf(fileHandle_p, " %f", sampleAvg[count]);          fprintf(fileHandle_p, " %e", sampleAvg[count]);
533          count++;          count++;
         elems++;  
534            }            }
535            for (int n=0; n<3-shape; n++) {            for (int n=0; n<3-shape; n++)
536          fprintf(fileHandle_p, " 0");              fprintf(fileHandle_p, " %e", 0.);
         elems++;  
           }  
         }  
         for (int m=0; m<nCompReqd-elems; m++) {  
           fprintf(fileHandle_p, " 0");  
537          }          }
538            for (int m=0; m<3-shape; m++) {
539                for (int n=0; n<3; n++)
540                   fprintf(fileHandle_p, " %e", 0.);
541                }
542        }        }
543        fprintf(fileHandle_p, "\n");        fprintf(fileHandle_p, "\n");
544      }      }
# Line 543  void Finley_Mesh_saveVTK(const char * fi Line 554  void Finley_Mesh_saveVTK(const char * fi
554            "NumberOfComponents=\"%d\" "            "NumberOfComponents=\"%d\" "
555            "format=\"ascii\">\n",            "format=\"ascii\">\n",
556            dataNameStr, nCompReqd);            dataNameStr, nCompReqd);
557        /* write out the data */
558        /* if the number of required components is more than the number
559         * of actual components, pad with zeros
560           */
561          bool_t do_write=TRUE;
562          int shape = getDataPointShape(data_p, 0);
563          if (shape > 3) {
564          sprintf(error_msg,"shape should be 1, 2, or 3; I got %d\n", shape);
565              Finley_setError(VALUE_ERROR,error_msg);
566          return;
567          }
568        for (i=0; i<mesh_p->Nodes->numNodes; i++) {        for (i=0; i<mesh_p->Nodes->numNodes; i++) {
569      switch (nodetype) {      switch (nodetype) {
570      case(FINLEY_DEGREES_OF_FREEDOM):      case(FINLEY_DEGREES_OF_FREEDOM):
571        values = getSampleData(data_p,        values = getSampleData(data_p,mesh_p->Nodes->degreeOfFreedom[i]);
                  mesh_p->Nodes->degreeOfFreedom[i]);  
572        break;        break;
573      case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):      case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):
574        if (mesh_p->Nodes->toReduced[i]>=0) {        if (mesh_p->Nodes->toReduced[i]>=0) {
575          values = getSampleData(data_p,              do_write=TRUE;
576                     mesh_p->Nodes->reducedDegreeOfFreedom[i]);          values = getSampleData(data_p,mesh_p->Nodes->reducedDegreeOfFreedom[i]);
577        }        } else {
578                do_write=FALSE;
579              }
580        break;        break;
581      case(FINLEY_NODES):      case(FINLEY_NODES):
582        values = getSampleData(data_p,i);        values = getSampleData(data_p,i);
583        break;        break;
584      }      }
     /* 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;  
     }  
585      /* write the data different ways for scalar, vector and tensor */      /* write the data different ways for scalar, vector and tensor */
586      if (nCompReqd == 1) {          if (do_write) {
587        fprintf(fileHandle_p, " %f", values[0]);          if (nCompReqd == 1) {
588      }            fprintf(fileHandle_p, " %e", 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++;  
589          }          }
590          for (int n=0; n<3-shape; n++) {          else if (nCompReqd == 3) {
591            fprintf(fileHandle_p, " 0");            /* write out the data */
592            elems++;            for (int m=0; m<shape; m++)
593                fprintf(fileHandle_p, " %e", values[m]);
594              /* pad with zeros */
595              for (int m=0; m<nCompReqd-shape; m++)
596                  fprintf(fileHandle_p, " %e", 0.);
597          }          }
598        }          else if (nCompReqd == 9) {
599        for (int m=0; m<nCompReqd-elems; m++) {            /* tensor data, so have a 3x3 matrix to output as a row
600          fprintf(fileHandle_p, " 0");             * of 9 data points */
601        }            int count = 0;
602      }            for (int m=0; m<shape; m++) {
603      fprintf(fileHandle_p, "\n");              for (int n=0; n<shape; n++) {
604                  fprintf(fileHandle_p, " %e", values[count]);
605                  count++;
606                }
607                for (int n=0; n<3-shape; n++)
608                  fprintf(fileHandle_p, " %e", 0.);
609              }
610              for (int m=0; m<3-shape; m++)  {
611                  for (int n=0; n<3; n++)
612                      fprintf(fileHandle_p, " %e", 0.);
613                  }
614            }
615            fprintf(fileHandle_p, "\n");
616             }
617        }        }
618        fprintf(fileHandle_p, "</DataArray>\n");        fprintf(fileHandle_p, "</DataArray>\n");
619        fprintf(fileHandle_p, "</PointData>\n");        fprintf(fileHandle_p, "</PointData>\n");
# Line 625  void Finley_Mesh_saveVTK(const char * fi Line 632  void Finley_Mesh_saveVTK(const char * fi
632  }  }
633    
634  /*  /*
  * $Log$  
635   * Revision 1.6  2005/08/12 01:45:43  jgs   * Revision 1.6  2005/08/12 01:45:43  jgs
636   * erge of development branch dev-02 back to main trunk on 2005-08-12   * erge of development branch dev-02 back to main trunk on 2005-08-12
637   *   *
638     * Revision 1.5.2.4  2005/09/09 08:15:17  gross
639     * bugs in vtk and dx writer fixed
640     *
641     * Revision 1.5.2.3  2005/09/08 08:28:39  gross
642     * some cleanup in savevtk
643     *
644     * Revision 1.5.2.2  2005/09/07 06:26:20  gross
645     * the solver from finley are put into the standalone package paso now
646     *
647   * Revision 1.5.2.1  2005/08/10 06:14:37  gross   * Revision 1.5.2.1  2005/08/10 06:14:37  gross
648   * QUADRATIC HEXAHEDRON elements fixed   * QUADRATIC HEXAHEDRON elements fixed
649   *   *

Legend:
Removed from v.147  
changed lines
  Added in v.150

  ViewVC Help
Powered by ViewVC 1.1.26