/[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

trunk/esys2/finley/src/finleyC/Mesh_saveVTK.c revision 147 by jgs, Fri Aug 12 01:45:47 2005 UTC trunk/finley/src/Mesh_saveVTK.c revision 1705 by ksteube, Thu Aug 14 05:56:40 2008 UTC
# Line 1  Line 1 
1    
2  /* $Id$ */  /* $Id$ */
3    
4    /*******************************************************
5     *
6     *           Copyright 2003-2007 by ACceSS MNRF
7     *       Copyright 2007 by University of Queensland
8     *
9     *                http://esscc.uq.edu.au
10     *        Primary Business: Queensland, Australia
11     *  Licensed under the Open Software License version 3.0
12     *     http://www.opensource.org/licenses/osl-3.0.php
13     *
14     *******************************************************/
15    
16  /**************************************************************/  /**************************************************************/
17    
18  /*   writes data and mesh in a vtk file */  /*   writes data and mesh in a vtk file */
19    /*   nodal data needs to be given on FINLEY_NODES or FINLEY_REDUCED_NODES */
20    
21  /**************************************************************/  /**************************************************************/
22    
 /*   Copyrights by ACcESS, Australia 2004 */  
 /*   Author: Paul Cochrane, cochrane@esscc.uq.edu.au */  
23    
 /**************************************************************/  
   
 #include "Finley.h"  
 #include "Common.h"  
24  #include "Mesh.h"  #include "Mesh.h"
25  #include "escript/Data/DataC.h"  #include "Assemble.h"
26  #include "vtkCellType.h"  /* copied from vtk source directory !!! */  #include "vtkCellType.h"  /* copied from vtk source directory !!! */
27    #include "paso/PasoUtil.h"
28    
29  void Finley_Mesh_saveVTK(const char * filename_p, Finley_Mesh *mesh_p, escriptDataC* data_p) {  #define LEN_PRINTED_INT_FORMAT (9+1)
30    /* if there is no mesh we just return */  #define INT_FORMAT "%d "
31    if (mesh_p==NULL) return;  #define INT_NEWLINE_FORMAT "%d\n"
32    Finley_ElementFile* elements=NULL;  #define FLOAT_SCALAR_FORMAT "%12.6e\n"
33    #define FLOAT_VECTOR_FORMAT "%12.6e %12.6e %12.6e\n"
34    #define FLOAT_TENSOR_FORMAT "%12.6e %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e\n"
35    #define LEN_PRINTED_FLOAT_SCALAR_FORMAT (12+1)
36    #define LEN_PRINTED_FLOAT_VECTOR_FORMAT (3*(12+1)+1)
37    #define LEN_PRINTED_FLOAT_TENSOR_FORMAT (9*(12+1)+1)
38    #define NEWLINE "\n"
39    #define LEN_TMP_BUFFER LEN_PRINTED_FLOAT_TENSOR_FORMAT+(MAX_numNodes*LEN_PRINTED_INT_FORMAT+1)+2
40    #define NCOMP_MAX 9
41    #define __STRCAT(dest,chunk,dest_in_use)  \
42    {                  \
43      strcpy(&dest[dest_in_use], chunk); \
44      dest_in_use+=strlen(chunk); \
45    }
46    
47    void Finley_Mesh_saveVTK(const char * filename_p,
48                             Finley_Mesh *mesh_p,
49                             const dim_t num_data,
50                             char* *names_p,
51                             escriptDataC* *data_pp)
52    {
53    #ifdef USE_VTK
54      char error_msg[LenErrorMsg_MAX], *txt_buffer=NULL, tmp_buffer[LEN_TMP_BUFFER];
55      double sampleAvg[NCOMP_MAX], *values, rtmp;
56      size_t txt_buffer_in_use;
57      dim_t len_txt_buffer,  max_len_names;
58      FILE * fileHandle_p = NULL;
59      int mpi_size, i, j, cellType;
60      dim_t i_data;
61      dim_t nDim, globalNumPoints, numCells, globalNumCells, numVTKNodesPerElement;
62      dim_t myNumPoints, numPointsPerSample, rank, nComp, nCompReqd;
63      dim_t shape, NN, numCellFactor, myNumCells, max_name_len;
64      bool_t *isCellCentered=NULL,write_celldata=FALSE,write_pointdata=FALSE;
65      bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
66      index_t myFirstNode, myLastNode, *globalNodeIndex, k, *node_index, myFirstCell;
67      #ifdef PASO_MPI
68      int ierr;
69      /* int amode = MPI_MODE_CREATE | MPI_MODE_WRONLY |  MPI_MODE_SEQUENTIAL;  */
70      const int amode = MPI_MODE_CREATE | MPI_MODE_WRONLY | MPI_MODE_UNIQUE_OPEN;
71      MPI_File mpi_fileHandle_p;
72      MPI_Status mpi_status;
73      MPI_Request mpi_req;
74      MPI_Info mpi_info=MPI_INFO_NULL;
75      #endif
76      Paso_MPI_rank my_mpi_rank;
77      int nodetype=FINLEY_NODES;
78      int elementtype=FINLEY_UNKNOWN;
79    char elemTypeStr[32];    char elemTypeStr[32];
80    int i, j, k, numVTKNodesPerElement, isCellCentered=FALSE, nodetype=FINLEY_DEGREES_OF_FREEDOM;    Finley_NodeMapping *nodeMapping=NULL;
81    double* values, rtmp;    Finley_ElementFile* elements=NULL;
82    int nDim = mesh_p->Nodes->numDim;    ElementTypeId TypeId;
83      
84    /* get a pointer to the relevant mesh component */  
85    if (isEmpty(data_p)) {    /****************************************/
86      elements=mesh_p->Elements;    /*                                      */
87    } else {    /*       tags in the vtk file           */
88      switch(getFunctionSpaceType(data_p)) {  
89      case(FINLEY_DEGREES_OF_FREEDOM):    char* tags_header="<?xml version=\"1.0\"?>\n" \
90        nodetype = FINLEY_DEGREES_OF_FREEDOM;                      "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\">\n" \
91        isCellCentered = FALSE;                      "<UnstructuredGrid>\n" \
92        elements = mesh_p->Elements;                      "<Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n" \
93        break;                      "<Points>\n" \
94      case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):                      "<DataArray NumberOfComponents=\"%d\" type=\"Float64\" format=\"ascii\">\n";
95        Finley_ErrorCode=VALUE_ERROR;    char *tag_End_DataArray = "</DataArray>\n";
96        sprintf(Finley_ErrorMsg,    char* tag_End_PointData =  "</PointData>\n";
97            "Reduced degrees of freedom is not yet "    char* tag_End_CellData =  "</CellData>\n";
98            "implemented for saving vtk files\n");    char *footer = "</Piece>\n</UnstructuredGrid>\n</VTKFile>\n";
99        return;    char* tags_End_Points_and_Start_Conn = "</DataArray>\n</Points>\n<Cells>\n<DataArray Name=\"connectivity\" type=\"Int32\" format=\"ascii\">\n" ;
100      case(FINLEY_NODES):    char* tags_End_Conn_and_Start_Offset = "</DataArray>\n<DataArray Name=\"offsets\" type=\"Int32\" format=\"ascii\">\n";
101        nodetype=FINLEY_NODES;    char* tags_End_Offset_and_Start_Type = "</DataArray>\n<DataArray Name=\"types\" type=\"UInt8\" format=\"ascii\">\n";
102        isCellCentered=FALSE;    char* tag_Float_DataArray="<DataArray Name=\"%s\" type=\"Float64\" NumberOfComponents=\"%d\" format=\"ascii\">\n";
103        elements=mesh_p->Elements;    char *tags_End_Type_And_Cells = "</DataArray>\n</Cells>\n";
       break;  
     case(FINLEY_ELEMENTS):  
       isCellCentered=TRUE;  
       elements=mesh_p->Elements;  
       break;  
     case(FINLEY_FACE_ELEMENTS):  
       isCellCentered=TRUE;  
       elements=mesh_p->FaceElements;  
       break;  
     case(FINLEY_POINTS):  
       isCellCentered=TRUE;  
       elements=mesh_p->Points;  
       break;  
     case(FINLEY_CONTACT_ELEMENTS_1):  
     case(FINLEY_CONTACT_ELEMENTS_2):  
       isCellCentered=TRUE;  
       elements=mesh_p->ContactElements;  
       break;  
     default:  
       Finley_ErrorCode=TYPE_ERROR;  
       sprintf(Finley_ErrorMsg,  
           "Finley does not know anything about function space type %d",  
           getFunctionSpaceType(data_p));  
       return;  
     }  
   }  
104    
105    /* the number of points */    int VTK_QUADRATIC_HEXAHEDRON_INDEX[] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 16, 17, 18, 19, 12, 13, 14, 15 };
106    int numPoints = mesh_p->Nodes->numNodes;    /* if there is no mesh we just return */
107      if (mesh_p==NULL) return;
108    
109    /* the number of cells */    my_mpi_rank = mesh_p->Nodes->MPIInfo->rank;
110    if (elements == NULL) {    mpi_size  = mesh_p->Nodes->MPIInfo->size;
111      Finley_ErrorCode = VALUE_ERROR;    nDim = mesh_p->Nodes->numDim;
112      sprintf(Finley_ErrorMsg,  
113          "elements object is NULL; cannot proceed");    if (! ( (nDim ==2) || (nDim == 3) ) ) {
114      return;          Finley_setError(IO_ERROR, "saveVTK: spatial dimension 2 or 3 is supported only.");
115            return;  
116    }    }
117    int numCells = elements->numElements;      /*************************************************************************************/
118      
119    /* open the file and check handle */    /* open the file and check handle */
120    FILE * fileHandle_p = fopen(filename_p, "w");  
121    if (fileHandle_p==NULL) {    if (mpi_size > 1) {
122      Finley_ErrorCode=IO_ERROR;          #ifdef PASO_MPI
123      sprintf(Finley_ErrorMsg,            /* Collective Call */
124          "File %s could not be opened for writing.", filename_p);            #ifdef MPIO_HINTS
125      return;              MPI_Info_create(&mpi_info);
126    }              /*  MPI_Info_set(mpi_info, "striping_unit",        "424288"); */
127    /* xml header */              /*  MPI_Info_set(mpi_info, "striping_factor",      "16"); */
128    fprintf(fileHandle_p, "<?xml version=\"1.0\"?>\n");              /*  MPI_Info_set(mpi_info, "collective_buffering", "true"); */
129    fprintf(fileHandle_p,              /*  MPI_Info_set(mpi_info, "cb_block_size",        "131072"); */
130        "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\">\n");              /*  MPI_Info_set(mpi_info, "cb_buffer_size",       "1048567"); */
131                /*  MPI_Info_set(mpi_info, "cb_nodes",             "8"); */
132    /* finley uses an unstructured mesh, so UnstructuredGrid *should* work */              /*    MPI_Info_set(mpi_info, "access_style", "write_once, sequential"); */
133    fprintf(fileHandle_p, "<UnstructuredGrid>\n");            
134                /*XFS only */
135    /* is there only one "piece" to the data?? */              /*   MPI_Info_set(mpi_info, "direct_write",          "true"); */
136    fprintf(fileHandle_p, "<Piece "            #endif
137        "NumberOfPoints=\"%d\" "            if ( my_mpi_rank == 0) {
138        "NumberOfCells=\"%d\">\n",                if  (Paso_fileExists(filename_p)) remove(filename_p);
139        numPoints, numCells);            }
140              ierr=MPI_File_open(mesh_p->Nodes->MPIInfo->comm, (char*)filename_p, amode,mpi_info, &mpi_fileHandle_p);
141    /* now for the points; equivalent to positions section in saveDX() */            if (ierr != MPI_SUCCESS) {
142    /* "The points element explicitly defines coordinates for each point            perror(filename_p);
143     * individually.  It contains one DataArray element describing an array                sprintf(error_msg, "saveVTK: File %s could not be opened for writing in parallel.", filename_p);
144     * with three components per value, each specifying the coordinates of one                Finley_setError(IO_ERROR,error_msg);
145     * point" - from Vtk User's Guide            } else {
146     */               MPI_File_set_view(mpi_fileHandle_p,MPI_DISPLACEMENT_CURRENT,MPI_CHAR, MPI_CHAR, "native" , mpi_info);
147    fprintf(fileHandle_p, "<Points>\n");            }
148    /*          #endif
    * 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");  
149    } else {    } else {
150      fprintf(fileHandle_p, "<DataArray "          fileHandle_p = fopen(filename_p, "w");
151          "NumberOfComponents=\"%d\" "          if (fileHandle_p==NULL) {
152          "type=\"Float32\" "             sprintf(error_msg, "saveVTK: File %s could not be opened for writing.", filename_p);
153          "format=\"ascii\">\n",             Finley_setError(IO_ERROR,error_msg);
154          nDim);           }
155    }    }
156    for (i = 0; i < mesh_p->Nodes->numNodes; i++) {    if (! Paso_MPIInfo_noError(mesh_p->Nodes->MPIInfo) ) return;
157      fprintf(fileHandle_p,    /*************************************************************************************/
158          "%e", mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)]);  
159      for (j = 1; j < nDim; j++) {    /* find the mesh type to be written */
       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");  
   
   /* 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 */  
   int cellType;  
   if (elements!=NULL) {  
     fprintf(fileHandle_p, "<Cells>\n");  
     ElementTypeId TypeId = elements->ReferenceElement->Type->TypeId;  
     switch(TypeId) {  
     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;  
     case Line2Face:  
       cellType = VTK_VERTEX;  
       break;  
     case Line3Face:  
       cellType = VTK_VERTEX;  
       break;  
     case Tri3Face:  
       cellType = VTK_LINE;  
       break;  
     case Tri6Face:  
       cellType = VTK_QUADRATIC_EDGE;  
       break;  
     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;  
     case Line2_Contact:  
       cellType = VTK_LINE;  
       break;  
     case Line3_Contact:  
       cellType = VTK_QUADRATIC_EDGE;  
       break;  
     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;  
     case Tri3Face_Contact:  
       cellType = VTK_LINE;  
       break;  
     case Tri6Face_Contact:  
       cellType = VTK_QUADRATIC_EDGE;  
       break;  
     case Rec4Face_Contact:  
       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:  
       numVTKNodesPerElement = 2;  
       strcpy(elemTypeStr, "VTK_LINE");  
       break;  
     case VTK_TRIANGLE:  
       numVTKNodesPerElement = 3;  
       strcpy(elemTypeStr, "VTK_TRIANGLE");  
       break;  
     case VTK_QUAD:  
       numVTKNodesPerElement = 4;  
       strcpy(elemTypeStr, "VTK_QUAD");  
       break;  
     case VTK_TETRA:  
       numVTKNodesPerElement = 4;  
       strcpy(elemTypeStr, "VTK_TETRA");  
       break;  
     case VTK_HEXAHEDRON:  
       numVTKNodesPerElement = 8;  
       strcpy(elemTypeStr, "VTK_HEXAHEDRON");  
       break;  
     case VTK_QUADRATIC_EDGE:  
       numVTKNodesPerElement = 3;  
       strcpy(elemTypeStr, "VTK_QUADRATIC_EDGE");  
       break;  
     case VTK_QUADRATIC_TRIANGLE:  
       numVTKNodesPerElement = 6;  
       strcpy(elemTypeStr, "VTK_QUADRATIC_TRIANGLE");  
       break;  
     case VTK_QUADRATIC_QUAD:  
       numVTKNodesPerElement = 8;  
       strcpy(elemTypeStr, "VTK_QUADRATIC_QUAD");  
       break;  
     case VTK_QUADRATIC_TETRA:  
       numVTKNodesPerElement = 10;  
       strcpy(elemTypeStr, "VTK_QUADRATIC_TETRA");  
       break;  
     case VTK_QUADRATIC_HEXAHEDRON:  
       numVTKNodesPerElement = 20;  
       strcpy(elemTypeStr, "VTK_QUADRATIC_HEXAHEDRON");  
       break;  
     default:  
       Finley_ErrorCode = VALUE_ERROR;  
       sprintf(Finley_ErrorMsg,  
           "Cell type %d is not supported by VTK", cellType);  
       return;  
     }  
160    
161      /* write out the DataArray element for the connectivity */    isCellCentered=TMPMEMALLOC(num_data,bool_t);
162      fprintf(fileHandle_p, "<DataArray "    max_len_names=0;
163          "Name=\"connectivity\" "    if (!Finley_checkPtr(isCellCentered)) {
164          "type=\"Int32\" "       nodetype=FINLEY_UNKNOWN;
165          "format=\"ascii\">\n");       elementtype=FINLEY_UNKNOWN;
166      int NN = elements->ReferenceElement->Type->numNodes;       for (i_data=0;i_data<num_data;++i_data) {
167      if (VTK_QUADRATIC_HEXAHEDRON==cellType) {         if (! isEmpty(data_pp[i_data])) {
168         for (i = 0; i < numCells; i++) {           switch(getFunctionSpaceType(data_pp[i_data]) ) {
169           fprintf(fileHandle_p,"%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d",           case FINLEY_NODES:
170                              mesh_p->Elements->Nodes[INDEX2(0, i, NN)],             nodetype = (nodetype == FINLEY_REDUCED_NODES) ? FINLEY_REDUCED_NODES : FINLEY_NODES;
171                              mesh_p->Elements->Nodes[INDEX2(1, i, NN)],             isCellCentered[i_data]=FALSE;
172                              mesh_p->Elements->Nodes[INDEX2(2, i, NN)],             if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
173                              mesh_p->Elements->Nodes[INDEX2(3, i, NN)],               elementtype=FINLEY_ELEMENTS;
174                              mesh_p->Elements->Nodes[INDEX2(4, i, NN)],             } else {
175                              mesh_p->Elements->Nodes[INDEX2(5, i, NN)],               Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
176                              mesh_p->Elements->Nodes[INDEX2(6, i, NN)],             }
177                              mesh_p->Elements->Nodes[INDEX2(7, i, NN)],             break;
178                              mesh_p->Elements->Nodes[INDEX2(8, i, NN)],           case FINLEY_REDUCED_NODES:
179                              mesh_p->Elements->Nodes[INDEX2(9, i, NN)],             nodetype = FINLEY_REDUCED_NODES;
180                              mesh_p->Elements->Nodes[INDEX2(10, i, NN)],             isCellCentered[i_data]=FALSE;
181                              mesh_p->Elements->Nodes[INDEX2(11, i, NN)],             if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
182                              mesh_p->Elements->Nodes[INDEX2(16, i, NN)],               elementtype=FINLEY_ELEMENTS;
183                              mesh_p->Elements->Nodes[INDEX2(17, i, NN)],             } else {
184                              mesh_p->Elements->Nodes[INDEX2(18, i, NN)],               Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
185                              mesh_p->Elements->Nodes[INDEX2(19, i, NN)],             }
186                              mesh_p->Elements->Nodes[INDEX2(12, i, NN)],             break;
187                              mesh_p->Elements->Nodes[INDEX2(13, i, NN)],           case FINLEY_ELEMENTS:
188                              mesh_p->Elements->Nodes[INDEX2(14, i, NN)],           case FINLEY_REDUCED_ELEMENTS:
189                              mesh_p->Elements->Nodes[INDEX2(15, i, NN)]);             isCellCentered[i_data]=TRUE;
190           fprintf(fileHandle_p, "\n");             if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
191                 elementtype=FINLEY_ELEMENTS;
192               } else {
193                 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
194               }
195               break;
196             case FINLEY_FACE_ELEMENTS:
197             case FINLEY_REDUCED_FACE_ELEMENTS:
198               isCellCentered[i_data]=TRUE;
199               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_FACE_ELEMENTS) {
200                 elementtype=FINLEY_FACE_ELEMENTS;
201               } else {
202                 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
203               }
204               break;
205             case FINLEY_POINTS:
206               isCellCentered[i_data]=TRUE;
207               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_POINTS) {
208                 elementtype=FINLEY_POINTS;
209               } else {
210                 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
211               }
212               break;
213             case FINLEY_CONTACT_ELEMENTS_1:
214             case FINLEY_REDUCED_CONTACT_ELEMENTS_1:
215               isCellCentered[i_data]=TRUE;
216               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {
217                 elementtype=FINLEY_CONTACT_ELEMENTS_1;
218               } else {
219                 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
220               }
221               break;
222             case FINLEY_CONTACT_ELEMENTS_2:
223             case FINLEY_REDUCED_CONTACT_ELEMENTS_2:
224               isCellCentered[i_data]=TRUE;
225               if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {
226                 elementtype=FINLEY_CONTACT_ELEMENTS_1;
227               } else {
228                 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
229               }
230               break;
231             default:
232               sprintf(error_msg,"saveVTK: unknown function space type %d",getFunctionSpaceType(data_pp[i_data]));
233               Finley_setError(TYPE_ERROR,error_msg);
234             }
235             if (isCellCentered[i_data]) {
236               write_celldata=TRUE;
237             } else {
238               write_pointdata=TRUE;
239             }
240             max_len_names =MAX(max_len_names,(dim_t)strlen(names_p[i_data]));
241         }         }
242      } else {       }
243         for (i = 0; i < numCells; i++) {       nodetype = (nodetype == FINLEY_UNKNOWN) ? FINLEY_NODES : nodetype;
244           for (j = 0; j < numVTKNodesPerElement; j++) fprintf(fileHandle_p,"%d ", mesh_p->Elements->Nodes[INDEX2(j, i, NN)]);    }
245           fprintf(fileHandle_p, "\n");    if (Finley_noError()) {
246    
247         /***************************************/
248    
249         /* select number of points and the mesh component */
250      
251         if (nodetype == FINLEY_REDUCED_NODES) {
252            myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
253            myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
254            globalNumPoints = Finley_NodeFile_getGlobalNumReducedNodes(mesh_p->Nodes);
255            globalNodeIndex= Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
256         } else {
257            myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
258            myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
259            globalNumPoints = Finley_NodeFile_getGlobalNumNodes(mesh_p->Nodes);
260            globalNodeIndex= Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
261         }
262         myNumPoints = myLastNode - myFirstNode;
263         if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;
264         switch(elementtype) {
265           case FINLEY_ELEMENTS:
266              elements=mesh_p->Elements;
267              break;
268            case FINLEY_FACE_ELEMENTS:
269              elements=mesh_p->FaceElements;
270              break;
271            case FINLEY_POINTS:
272              elements=mesh_p->Points;
273              break;
274            case FINLEY_CONTACT_ELEMENTS_1:
275              elements=mesh_p->ContactElements;
276              break;
277         }
278         if (elements==NULL) {
279           Finley_setError(SYSTEM_ERROR,"saveVTK: undefined element file");
280         } else {
281           /* map finley element type to VTK element type */
282           numCells = elements->numElements;
283           globalNumCells = Finley_ElementFile_getGlobalNumElements(elements);
284           myNumCells= Finley_ElementFile_getMyNumElements(elements);
285           myFirstCell= Finley_ElementFile_getFirstElement(elements);
286           NN = elements->numNodes;
287           if (nodetype==FINLEY_REDUCED_NODES) {
288              TypeId = elements->LinearReferenceElement->Type->TypeId;
289           } else {
290              TypeId = elements->ReferenceElement->Type->TypeId;
291         }         }
292      }         switch(TypeId) {
293      fprintf(fileHandle_p, "\n");          case Point1:
294      fprintf(fileHandle_p, "</DataArray>\n");          case Line2Face:
295            case Line3Face:
296      /* write out the DataArray element for the offsets */          case Point1_Contact:
297      fprintf(fileHandle_p, "<DataArray "          case Line2Face_Contact:
298          "Name=\"offsets\" "          case Line3Face_Contact:
299          "type=\"Int32\" "            numCellFactor=1;
300          "format=\"ascii\">\n");            cellType = VTK_VERTEX;
301      int counter = 0;  /* counter for the number of offsets written to file */            numVTKNodesPerElement = 1;
302      for (i=numVTKNodesPerElement; i<=numCells*numVTKNodesPerElement; i+=numVTKNodesPerElement) {            strcpy(elemTypeStr, "VTK_VERTEX");
303        fprintf(fileHandle_p, "%d ", i);            break;
304        counter++;        
305        /* if the counter gets too big (i.e. the line gets too long),          case Line2:
306         * then add a newline and reset */          case Tri3Face:
307        if (counter > 19) {          case Rec4Face:
308        counter = 0;          case Line2_Contact:
309        fprintf(fileHandle_p, "\n");          case Tri3_Contact:
310            case Tri3Face_Contact:
311            case Rec4Face_Contact:
312              numCellFactor=1;
313              cellType = VTK_LINE;
314              numVTKNodesPerElement = 2;
315              strcpy(elemTypeStr, "VTK_LINE");
316              break;
317          
318            case Tri3:
319            case Tet4Face:
320            case Tet4Face_Contact:
321              numCellFactor=1;
322              cellType = VTK_TRIANGLE;
323              numVTKNodesPerElement = 3;
324              strcpy(elemTypeStr, "VTK_TRIANGLE");
325              break;
326          
327            case Rec4:
328            case Hex8Face:
329            case Rec4_Contact:
330            case Hex8Face_Contact:
331              numCellFactor=1;
332              cellType = VTK_QUAD;
333              numVTKNodesPerElement = 4;
334              strcpy(elemTypeStr, "VTK_QUAD");
335              break;
336          
337            case Tet4:
338              numCellFactor=1;
339              cellType = VTK_TETRA;
340              numVTKNodesPerElement = 4;
341              strcpy(elemTypeStr, "VTK_TETRA");
342              break;
343          
344            case Hex8:
345              numCellFactor=1;
346              cellType = VTK_HEXAHEDRON;
347              numVTKNodesPerElement = 8;
348              strcpy(elemTypeStr, "VTK_HEXAHEDRON");
349              break;
350          
351            case Line3:
352            case Tri6Face:
353            case Rec8Face:
354            case Line3_Contact:
355            case Tri6Face_Contact:
356            case Rec8Face_Contact:
357              numCellFactor=1;
358              cellType = VTK_QUADRATIC_EDGE;
359              numVTKNodesPerElement = 3;
360              strcpy(elemTypeStr, "VTK_QUADRATIC_EDGE");
361              break;
362          
363            case Tri6:
364            case Tet10Face:
365            case Tri6_Contact:
366            case Tet10Face_Contact:
367              numCellFactor=1;
368              cellType = VTK_QUADRATIC_TRIANGLE;
369              numVTKNodesPerElement = 6;
370              strcpy(elemTypeStr, "VTK_QUADRATIC_TRIANGLE");
371              break;
372          
373            case Rec8:
374            case Hex20Face:
375            case Rec8_Contact:
376            case Hex20Face_Contact:
377              numCellFactor=1;
378              cellType = VTK_QUADRATIC_QUAD;
379              numVTKNodesPerElement = 8;
380              strcpy(elemTypeStr, "VTK_QUADRATIC_QUAD");
381              break;
382          
383            case Tet10:
384              numCellFactor=1;
385              cellType = VTK_QUADRATIC_TETRA;
386              numVTKNodesPerElement = 10;
387              strcpy(elemTypeStr, "VTK_QUADRATIC_TETRA");
388              break;
389          
390            case Hex20:
391              numCellFactor=1;
392              cellType = VTK_QUADRATIC_HEXAHEDRON;
393              numVTKNodesPerElement = 20;
394              strcpy(elemTypeStr, "VTK_QUADRATIC_HEXAHEDRON");
395              break;
396          
397            default:
398              sprintf(error_msg, "saveVTK: Element type %s is not supported by VTK",elements->ReferenceElement->Type->Name);
399              Finley_setError(VALUE_ERROR,error_msg);
400            }
401         }
402      }
403      /***************************************/
404    
405      /***************************************/
406      /*                                     */
407      /*   allocate text buffer              */
408      /*                                     */
409      max_name_len=0;
410      for (i_data =0 ;i_data<num_data;++i_data) max_name_len=MAX(max_name_len,(dim_t)strlen(names_p[i_data]));
411      len_txt_buffer= strlen(tags_header) + 3 * LEN_PRINTED_INT_FORMAT + (30+3*max_name_len); /* header */
412      if (mpi_size > 1) len_txt_buffer=MAX(len_txt_buffer, myNumPoints * LEN_TMP_BUFFER);
413      if (mpi_size > 1) len_txt_buffer=MAX(len_txt_buffer, numCellFactor*myNumCells*(LEN_PRINTED_INT_FORMAT*numVTKNodesPerElement+1));
414      len_txt_buffer=MAX(len_txt_buffer,200+3*max_len_names);
415      len_txt_buffer=MAX(len_txt_buffer, (dim_t)strlen(tag_Float_DataArray) + LEN_PRINTED_INT_FORMAT + max_len_names);
416      if (mpi_size > 1) len_txt_buffer=MAX(len_txt_buffer, numCellFactor*myNumCells*LEN_PRINTED_FLOAT_TENSOR_FORMAT);
417      if (mpi_size > 1) len_txt_buffer=MAX(len_txt_buffer, myNumPoints*LEN_PRINTED_FLOAT_TENSOR_FORMAT);
418      txt_buffer=TMPMEMALLOC(len_txt_buffer+1,char);
419      Finley_checkPtr(txt_buffer);
420      
421      if (Finley_noError()) {
422    
423         /* select number of points and the mesh component */
424    
425         sprintf(txt_buffer,tags_header,globalNumPoints,numCellFactor*globalNumCells,3);
426    
427          if (mpi_size > 1) {
428              if ( my_mpi_rank == 0) {
429                #ifdef PASO_MPI
430                  MPI_File_iwrite_shared(mpi_fileHandle_p,txt_buffer,strlen(txt_buffer),MPI_CHAR,&mpi_req);
431                  MPI_Wait(&mpi_req,&mpi_status);
432                #endif
433              }
434          } else {
435             fprintf(fileHandle_p,txt_buffer);
436        }        }
     }  
     fprintf(fileHandle_p, "\n");  
     fprintf(fileHandle_p, "</DataArray>\n");  
437    
438      /* write out the DataArray element for the types */        /* write the nodes */
439      counter = 0; /* counter for the number of types written to file */        
440      fprintf(fileHandle_p, "<DataArray "        if (mpi_size > 1) {
441          "Name=\"types\" "  
442          "type=\"UInt8\" "           txt_buffer[0] = '\0';
443          "format=\"ascii\">\n");           txt_buffer_in_use=0;
444      for (i=0; i<numCells; i++) {           if (nDim==2) {
445        fprintf(fileHandle_p, "%d ", cellType);              for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
446        counter++;                 if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
447        /* if the counter gets too big (i.e. the line gets too long),                   sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,
448         * then add a newline and reset */                                      mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
449        if (counter > 30) {                                      mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
450        counter = 0;                                      0.);
451        fprintf(fileHandle_p, "\n");                   __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use);
452                   }
453                }      
454             } else {
455                for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
456                   if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
457                     sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,
458                                                      mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
459                                                      mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
460                                                      mesh_p->Nodes->Coordinates[INDEX2(2, i, nDim)]);
461                     __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use);
462                   }
463                }    
464      
465             }
466             #ifdef PASO_MPI
467                MPI_File_write_ordered(mpi_fileHandle_p, txt_buffer,txt_buffer_in_use, MPI_CHAR, &mpi_status);
468             #endif    
469          } else {
470             if (nDim==2) {
471                for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
472                   if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
473                     fprintf(fileHandle_p,FLOAT_VECTOR_FORMAT,
474                                          mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
475                                          mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
476                                          0.);
477                   }
478                }      
479             } else {
480                for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
481                   if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
482                     fprintf(fileHandle_p,FLOAT_VECTOR_FORMAT,
483                                                  mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
484                                                  mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
485                                                  mesh_p->Nodes->Coordinates[INDEX2(2, i, nDim)]);
486                   }
487                }    
488      
489             }
490        }        }
     }  
     fprintf(fileHandle_p, "\n");  
     fprintf(fileHandle_p, "</DataArray>\n");  
491    
492      /* finish off the <Cells> element */        /* close the Points and open connectivity */
     fprintf(fileHandle_p, "</Cells>\n");  
   }  
493    
494    /* data */        if (mpi_size > 1) {
495    if (!isEmpty(data_p)) {            if ( my_mpi_rank == 0) {
496      int rank = getDataPointRank(data_p);               #ifdef PASO_MPI
497      int nComp = getDataPointSize(data_p);                  MPI_File_iwrite_shared(mpi_fileHandle_p, tags_End_Points_and_Start_Conn, strlen(tags_End_Points_and_Start_Conn), MPI_CHAR, &mpi_req);
498      /* barf if rank is greater than two */                  MPI_Wait(&mpi_req,&mpi_status);
499      if (rank > 2) {               #endif
500        Finley_ErrorCode = VALUE_ERROR;            }
501        sprintf(Finley_ErrorMsg,        } else {
502            "Vtk can't handle objects with rank greater than 2\n"           fprintf(fileHandle_p,tags_End_Points_and_Start_Conn);
503            "object rank = %d\n", rank);        }
504        return;  
505      }       /* write the cells */
506      /* if the rank == 0:   --> scalar data       if (nodetype == FINLEY_REDUCED_NODES) {
507       * if the rank == 1:   --> vector data          node_index=elements->ReferenceElement->Type->linearNodes;
508       * if the rank == 2:   --> tensor data       } else if (VTK_QUADRATIC_HEXAHEDRON==cellType) {
509       */          node_index=VTK_QUADRATIC_HEXAHEDRON_INDEX;
510      char dataNameStr[31], dataTypeStr[63];       } else if (numVTKNodesPerElement!=NN) {
511      int nCompReqd=1;   /* the number of components required by vtk */          node_index=elements->ReferenceElement->Type->geoNodes;
512      if (rank == 0) {       } else {
513        strcpy(dataNameStr, "escript_scalar_data");          node_index=NULL;
514        sprintf(dataTypeStr, "Scalars=\"%s\"", dataNameStr);       }
515        nCompReqd = 1;  
516      }       if ( mpi_size > 1) {
517      else if (rank == 1) {          txt_buffer[0] = '\0';
518        strcpy(dataNameStr, "escript_vector_data");          txt_buffer_in_use=0;
519        sprintf(dataTypeStr, "Vectors=\"%s\"", dataNameStr);          if (node_index == NULL) {
520        nCompReqd = 3;             for (i = 0; i < numCells; i++) {
521                  if (elements->Owner[i] == my_mpi_rank) {
522                     for (j = 0; j < numVTKNodesPerElement; j++) {
523                         sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(j, i, NN)]]);
524                         __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use)
525                     }
526                     __STRCAT(txt_buffer,NEWLINE,txt_buffer_in_use)
527                  }
528               }
529            } else {
530               for (i = 0; i < numCells; i++) {
531                  if (elements->Owner[i] == my_mpi_rank) {
532                     for (j = 0; j < numVTKNodesPerElement; j++) {
533                         sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(node_index[j], i, NN)]]);
534                         __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use)
535                     }
536                     __STRCAT(txt_buffer,NEWLINE,txt_buffer_in_use)
537                  }
538               }
539            }
540            #ifdef PASO_MPI
541               MPI_File_write_ordered(mpi_fileHandle_p,txt_buffer,txt_buffer_in_use, MPI_CHAR, &mpi_status);
542            #endif    
543         } else {
544            if (node_index == NULL) {
545               for (i = 0; i < numCells; i++) {
546                  for (j = 0; j < numVTKNodesPerElement; j++) {
547                     fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(j, i, NN)]]);
548                   }
549                  fprintf(fileHandle_p,NEWLINE);
550               }
551            } else {
552               for (i = 0; i < numCells; i++) {
553                  for (j = 0; j < numVTKNodesPerElement; j++) {
554                     fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(node_index[j], i, NN)]]);
555                   }
556                  fprintf(fileHandle_p,NEWLINE);
557               }
558            }
559    
560         }
561        
562         /* finalize the connection and start the offset section */
563         if (mpi_size > 1) {
564            if( my_mpi_rank == 0) {
565               #ifdef PASO_MPI
566                  MPI_File_iwrite_shared(mpi_fileHandle_p,tags_End_Conn_and_Start_Offset,strlen(tags_End_Conn_and_Start_Offset),MPI_CHAR,&mpi_req);
567                  MPI_Wait(&mpi_req,&mpi_status);
568               #endif
569            }
570         } else {
571            fprintf(fileHandle_p,tags_End_Conn_and_Start_Offset);
572         }
573    
574        /* write the offsets */
575          
576         if ( mpi_size > 1) {
577            txt_buffer[0] = '\0';
578            txt_buffer_in_use=0;
579            for (i=numVTKNodesPerElement*(myFirstCell*numCellFactor+1); i<=(myFirstCell+myNumCells)*numVTKNodesPerElement*numCellFactor; i+=numVTKNodesPerElement) {
580                sprintf(tmp_buffer, INT_NEWLINE_FORMAT, i);
581                __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use);
582             }
583             #ifdef PASO_MPI
584                MPI_File_write_ordered(mpi_fileHandle_p,txt_buffer,txt_buffer_in_use, MPI_CHAR, &mpi_status);
585             #endif    
586         } else {
587            for (i=numVTKNodesPerElement; i<=numCells*numVTKNodesPerElement*numCellFactor; i+=numVTKNodesPerElement) {
588               fprintf(fileHandle_p, INT_NEWLINE_FORMAT, i);
589            }
590        
591         }
592         /* finalize the offset section and start the type section */
593         if ( mpi_size > 1) {
594            if ( my_mpi_rank == 0) {
595               #ifdef PASO_MPI
596                  MPI_File_iwrite_shared(mpi_fileHandle_p,tags_End_Offset_and_Start_Type,strlen(tags_End_Offset_and_Start_Type),MPI_CHAR,&mpi_req);
597                  MPI_Wait(&mpi_req,&mpi_status);
598               #endif
599            }
600        } else {
601           fprintf(fileHandle_p,tags_End_Offset_and_Start_Type);
602      }      }
603      else if (rank == 2) {       /* write element type */
604        strcpy(dataNameStr, "escript_tensor_data");       sprintf(tmp_buffer, INT_NEWLINE_FORMAT, cellType);
605        sprintf(dataTypeStr, "Tensors=\"%s\"", dataNameStr);       if ( mpi_size > 1) {
606        nCompReqd = 9;          txt_buffer[0] = '\0';
607            txt_buffer_in_use=0;
608            for (i=0; i<numCells*numCellFactor; i++) __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use);
609             #ifdef PASO_MPI
610                MPI_File_write_ordered(mpi_fileHandle_p,txt_buffer,txt_buffer_in_use, MPI_CHAR, &mpi_status);
611             #endif    
612         } else {
613            for (i=0; i<numCells*numCellFactor; i++) fprintf(fileHandle_p, tmp_buffer);
614         }
615         /* finalize cell information */
616         if ( mpi_size > 1) {
617            if ( my_mpi_rank == 0) {
618               #ifdef PASO_MPI
619                  MPI_File_iwrite_shared(mpi_fileHandle_p,tags_End_Type_And_Cells,strlen(tags_End_Type_And_Cells),MPI_CHAR,&mpi_req);
620                  MPI_Wait(&mpi_req,&mpi_status);
621               #endif
622            }
623        } else {
624           fprintf(fileHandle_p,tags_End_Type_And_Cells);
625      }      }
626      /* if have cell centred data then use CellData element,   }
627       * if have node centred data, then use PointData element  
628       */   /* Write cell data */
629      if (isCellCentered) {   if (write_celldata && Finley_noError()) {
630        /* now for the cell data */        /* mark the active data arrays */
631        fprintf(fileHandle_p, "<CellData %s>\n", dataTypeStr);        txt_buffer[0] = '\0';
632        fprintf(fileHandle_p,        set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
633            "<DataArray "        strcat(txt_buffer, "<CellData");
634            "Name=\"%s\" "        for (i_data =0 ;i_data<num_data;++i_data) {
635            "type=\"Float32\" "          if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data]) {
636            "NumberOfComponents=\"%d\" "            /* if the rank == 0:   --> scalar data */
637            "format=\"ascii\">\n",            /* if the rank == 1:   --> vector data */
638            dataNameStr, nCompReqd);            /* if the rank == 2:   --> tensor data */
639        int numPointsPerSample = elements->ReferenceElement->numQuadNodes;  
640        if (numPointsPerSample) {            switch(getDataPointRank(data_pp[i_data])) {
641      for (i=0; i<numCells; i++) {            case 0:
642        values = getSampleData(data_p, i);              if (! set_scalar) {
643        double sampleAvg[nComp];                strcat(txt_buffer," Scalars=\"");
644        for (k=0; k<nComp; k++) {                strcat(txt_buffer,names_p[i_data]);
645          /* averaging over the number of points in the sample */                strcat(txt_buffer,"\"");
646          rtmp = 0.;                set_scalar=TRUE;
647          for (j=0; j<numPointsPerSample; j++) {              }
648            rtmp += values[INDEX2(k,j,nComp)];              break;
649          }            case 1:
650          sampleAvg[k] = rtmp/numPointsPerSample;              if (! set_vector) {
651        }                strcat(txt_buffer," Vectors=\"");
652        /* if the number of required components is more than the number                strcat(txt_buffer,names_p[i_data]);
653         * of actual components, pad with zeros                strcat(txt_buffer,"\"");
654         */                set_vector=TRUE;
655        /* probably only need to get shape of first element */              }
656        int shape = getDataPointShape(data_p, 0);              break;
657        if (shape > 3) {            case 2:
658          Finley_ErrorCode = VALUE_ERROR;              if (! set_tensor) {
659          sprintf(Finley_ErrorMsg,                strcat(txt_buffer," Tensors=\"");
660              "shape should be 1, 2, or 3; I got %d\n", shape);                strcat(txt_buffer,names_p[i_data]);
661          return;                strcat(txt_buffer,"\"");
662        }                set_tensor=TRUE;
663        /* write the data different ways for scalar, vector and tensor */              }
664        if (nCompReqd == 1) {              break;
665          fprintf(fileHandle_p, " %f", sampleAvg[0]);            default:
666        }              sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);
667        else if (nCompReqd == 3) {              Finley_setError(VALUE_ERROR,error_msg);
668          int shape = getDataPointShape(data_p, 0);              return;
669          /* write out the data */            }
670          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");  
     }  
671        }        }
672        fprintf(fileHandle_p, "</DataArray>\n");        strcat(txt_buffer, ">\n");
673        fprintf(fileHandle_p, "</CellData>\n");        if ( mpi_size > 1) {
674      } else {          if ( my_mpi_rank == 0) {
675        /* now for the point data */             #ifdef PASO_MPI
676        fprintf(fileHandle_p, "<PointData %s>\n", dataTypeStr);                MPI_File_iwrite_shared(mpi_fileHandle_p,txt_buffer,strlen(txt_buffer),MPI_CHAR,&mpi_req);
677        fprintf(fileHandle_p, "<DataArray "                MPI_Wait(&mpi_req,&mpi_status);
678            "Name=\"%s\" "             #endif
679            "type=\"Float32\" "          }
680            "NumberOfComponents=\"%d\" "        } else {
681            "format=\"ascii\">\n",            fprintf(fileHandle_p,txt_buffer);
682            dataNameStr, nCompReqd);        }
683        for (i=0; i<mesh_p->Nodes->numNodes; i++) {        /* write the arrays */
684      switch (nodetype) {        for (i_data =0 ;i_data<num_data;++i_data) {
685      case(FINLEY_DEGREES_OF_FREEDOM):           if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data]) {
686        values = getSampleData(data_p,              txt_buffer[0] = '\0';
687                   mesh_p->Nodes->degreeOfFreedom[i]);              txt_buffer_in_use=0;
688        break;              numPointsPerSample=getNumDataPointsPerSample(data_pp[i_data]);
689      case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):              rank = getDataPointRank(data_pp[i_data]);
690        if (mesh_p->Nodes->toReduced[i]>=0) {              nComp = getDataPointSize(data_pp[i_data]);
691          values = getSampleData(data_p,              nCompReqd=1;   /* the number of components mpi_required by vtk */
692                     mesh_p->Nodes->reducedDegreeOfFreedom[i]);              shape=0;
693        }              if (rank == 0) {
694        break;                nCompReqd = 1;
695      case(FINLEY_NODES):              } else if (rank == 1) {
696        values = getSampleData(data_p,i);                shape=getDataPointShape(data_pp[i_data], 0);
697        break;                if  (shape>3) {
698      }                  Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");
699      /* write out the data */                }
700      /* if the number of required components is more than the number                nCompReqd = 3;
701       * of actual components, pad with zeros              } else {
702       */                shape=getDataPointShape(data_pp[i_data], 0);
703      /* probably only need to get shape of first element */                if  (shape>3 || shape != getDataPointShape(data_pp[i_data], 1)) {
704      int shape = getDataPointShape(data_p, 0);                  Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
705      if (shape > 3) {                }
706        Finley_ErrorCode = VALUE_ERROR;                nCompReqd = 9;
707        sprintf(Finley_ErrorMsg,              }
708            "shape should be 1, 2, or 3; I got %d\n", shape);              if (Finley_noError()) {
709        return;                 sprintf(txt_buffer,tag_Float_DataArray,names_p[i_data], nCompReqd);
710      }                 if ( mpi_size > 1) {
711      /* write the data different ways for scalar, vector and tensor */                   if ( my_mpi_rank == 0) {
712      if (nCompReqd == 1) {                      #ifdef PASO_MPI
713        fprintf(fileHandle_p, " %f", values[0]);                         MPI_File_iwrite_shared(mpi_fileHandle_p,txt_buffer,strlen(txt_buffer),MPI_CHAR,&mpi_req);
714      }                         MPI_Wait(&mpi_req,&mpi_status);
715      else if (nCompReqd == 3) {                      #endif
716        int shape = getDataPointShape(data_p, 0);                   }
717        /* write out the data */                 } else {
718        for (int m=0; m<shape; m++) {                     fprintf(fileHandle_p,txt_buffer);
719          fprintf(fileHandle_p, " %f", values[m]);                 }
720        }                 for (i=0; i<numCells; i++) {
721        /* pad with zeros */                     if (elements->Owner[i] == my_mpi_rank) {
722        for (int m=0; m<nCompReqd-shape; m++) {                        values = getSampleData(data_pp[i_data], i);
723          fprintf(fileHandle_p, " 0");                        /* averaging over the number of points in the sample */
724        }                        for (k=0; k<MIN(nComp,NCOMP_MAX); k++) {
725      }                           if (isExpanded(data_pp[i_data])) {
726      else if (nCompReqd == 9) {                             rtmp = 0.;
727        /* tensor data, so have a 3x3 matrix to output as a row                             for (j=0; j<numPointsPerSample; j++) rtmp += values[INDEX2(k,j,nComp)];
728         * of 9 data points */                             sampleAvg[k] = rtmp/numPointsPerSample;
729        int count = 0;                          } else {
730        int elems = 0;                             sampleAvg[k] = values[k];
731        for (int m=0; m<shape; m++) {                          }
732          for (int n=0; n<shape; n++) {                        }
733            fprintf(fileHandle_p, " %f", values[count]);                        /* if the number of mpi_required components is more than the number
734            count++;                        * of actual components, pad with zeros
735            elems++;                        */
736          }                        /* probably only need to get shape of first element */
737          for (int n=0; n<3-shape; n++) {                        /* write the data different ways for scalar, vector and tensor */
738            fprintf(fileHandle_p, " 0");                        if (nCompReqd == 1) {
739            elems++;                          sprintf(tmp_buffer,FLOAT_SCALAR_FORMAT,sampleAvg[0]);
740          }                        } else if (nCompReqd == 3) {
741        }                          if (shape==1) {
742        for (int m=0; m<nCompReqd-elems; m++) {                           sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,sampleAvg[0],0.,0.);
743          fprintf(fileHandle_p, " 0");                          } else if (shape==2) {
744        }                           sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,sampleAvg[0],sampleAvg[1],0.);
745      }                          } else if (shape==3) {
746      fprintf(fileHandle_p, "\n");                           sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,sampleAvg[0],sampleAvg[1],sampleAvg[2]);
747                            }
748                          } else if (nCompReqd == 9) {
749                            if (shape==1) {
750                             sprintf(tmp_buffer,FLOAT_TENSOR_FORMAT,sampleAvg[0],0.,0.,
751                                                                    0.,0.,0.,
752                                                                    0.,0.,0.);
753                            } else if (shape==2) {
754                             sprintf(tmp_buffer,FLOAT_TENSOR_FORMAT,sampleAvg[0],sampleAvg[1],0.,
755                                                                    sampleAvg[2],sampleAvg[3],0.,
756                                                                    0.,0.,0.);
757                            } else if (shape==3) {
758                             sprintf(tmp_buffer,FLOAT_TENSOR_FORMAT,sampleAvg[0],sampleAvg[1],sampleAvg[2],
759                                                                    sampleAvg[3],sampleAvg[4],sampleAvg[5],
760                                                                    sampleAvg[6],sampleAvg[7],sampleAvg[8]);
761                            }
762                          }
763                          if ( mpi_size > 1) {
764                            __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use);
765                          } else {
766                            fprintf(fileHandle_p,tmp_buffer);
767                          }
768                      }
769                   }
770                   if ( mpi_size > 1) {
771                         #ifdef PASO_MPI
772                            MPI_File_write_ordered(mpi_fileHandle_p,txt_buffer,txt_buffer_in_use, MPI_CHAR, &mpi_status);
773                         #endif    
774                         if ( my_mpi_rank == 0) {
775                            #ifdef PASO_MPI
776                               MPI_File_iwrite_shared(mpi_fileHandle_p,tag_End_DataArray,strlen(tag_End_DataArray),MPI_CHAR,&mpi_req);
777                               MPI_Wait(&mpi_req,&mpi_status);
778                            #endif
779                         }
780                   } else {
781                       fprintf(fileHandle_p,tag_End_DataArray);
782                   }
783                }
784             }
785          }
786          if ( mpi_size > 1) {
787            if ( my_mpi_rank == 0) {
788               #ifdef PASO_MPI
789                  MPI_File_iwrite_shared(mpi_fileHandle_p,tag_End_CellData,strlen(tag_End_CellData),MPI_CHAR,&mpi_req);
790                  MPI_Wait(&mpi_req,&mpi_status);
791               #endif
792            }
793          } else {
794              fprintf(fileHandle_p,tag_End_CellData);
795        }        }
       fprintf(fileHandle_p, "</DataArray>\n");  
       fprintf(fileHandle_p, "</PointData>\n");  
     }  
796    }    }
797      /* point data */
798    /* finish off the piece */    if (write_pointdata && Finley_noError()) {
799    fprintf(fileHandle_p, "</Piece>\n");        /* mark the active data arrays */
800          set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
801    fprintf(fileHandle_p, "</UnstructuredGrid>\n");        txt_buffer[0] = '\0';
802    /* write the xml footer */        strcat(txt_buffer, "<PointData");
803    fprintf(fileHandle_p, "</VTKFile>\n");        for (i_data =0 ;i_data<num_data;++i_data) {
804    /* close the file */          if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data]) {
805    fclose(fileHandle_p);            /* if the rank == 0:   --> scalar data */
806              /* if the rank == 1:   --> vector data */
807              /* if the rank == 2:   --> tensor data */
808    
809              switch(getDataPointRank(data_pp[i_data])) {
810              case 0:
811                if (! set_scalar) {
812                  strcat(txt_buffer," Scalars=\"");
813                  strcat(txt_buffer,names_p[i_data]);
814                  strcat(txt_buffer,"\"");
815                  set_scalar=TRUE;
816                }
817                break;
818              case 1:
819                if (! set_vector) {
820                  strcat(txt_buffer," Vectors=\"");
821                  strcat(txt_buffer,names_p[i_data]);
822                  strcat(txt_buffer,"\"");
823                  set_vector=TRUE;
824                }
825                break;
826              case 2:
827                if (! set_tensor) {
828                  strcat(txt_buffer," Tensors=\"");
829                  strcat(txt_buffer,names_p[i_data]);
830                  strcat(txt_buffer,"\"");
831                  set_tensor=TRUE;
832                }
833                break;
834              default:
835                sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);
836                Finley_setError(VALUE_ERROR,error_msg);
837                return;
838              }
839            }
840          }
841          strcat(txt_buffer, ">\n");
842          if ( mpi_size > 1) {
843            if ( my_mpi_rank == 0) {
844               #ifdef PASO_MPI
845                  MPI_File_iwrite_shared(mpi_fileHandle_p,txt_buffer,strlen(txt_buffer),MPI_CHAR,&mpi_req);
846                  MPI_Wait(&mpi_req,&mpi_status);
847               #endif
848            }
849          } else {
850              fprintf(fileHandle_p,txt_buffer);
851          }
852          /* write the arrays */
853          for (i_data =0 ;i_data<num_data;++i_data) {
854             if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data]) {
855                txt_buffer[0] = '\0';
856                txt_buffer_in_use=0;
857                numPointsPerSample=getNumDataPointsPerSample(data_pp[i_data]);
858                rank = getDataPointRank(data_pp[i_data]);
859                nComp = getDataPointSize(data_pp[i_data]);
860                if (getFunctionSpaceType(data_pp[i_data]) == FINLEY_REDUCED_NODES) {
861                   nodeMapping=mesh_p->Nodes->reducedNodesMapping;
862                } else {
863                   nodeMapping=mesh_p->Nodes->nodesMapping;
864                }
865                nCompReqd=1;   /* the number of components mpi_required by vtk */
866                shape=0;
867                if (rank == 0) {
868                  nCompReqd = 1;
869                } else if (rank == 1) {
870                  shape=getDataPointShape(data_pp[i_data], 0);
871                  if  (shape>3) {
872                    Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");
873                  }
874                  nCompReqd = 3;
875                } else {
876                  shape=getDataPointShape(data_pp[i_data], 0);
877                  if  (shape>3 || shape != getDataPointShape(data_pp[i_data], 1)) {
878                    Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
879                  }
880                  nCompReqd = 9;
881                }
882                if (Finley_noError()) {
883                   sprintf(txt_buffer,tag_Float_DataArray,names_p[i_data], nCompReqd);
884                   if ( mpi_size > 1) {
885                     if ( my_mpi_rank == 0) {
886                        #ifdef PASO_MPI
887                           MPI_File_iwrite_shared(mpi_fileHandle_p,txt_buffer,strlen(txt_buffer),MPI_CHAR,&mpi_req);
888                           MPI_Wait(&mpi_req,&mpi_status);
889                        #endif
890                     }
891                   } else {
892                       fprintf(fileHandle_p,txt_buffer);
893                   }
894                   for (i=0; i<mesh_p->Nodes->numNodes; i++) {
895                      k=globalNodeIndex[i];
896                      if ( (myFirstNode <= k) && (k < myLastNode) ) {
897                         values = getSampleData(data_pp[i_data], nodeMapping->target[i]);
898                         /* if the number of mpi_required components is more than the number
899                         * of actual components, pad with zeros
900                         */
901                         /* probably only need to get shape of first element */
902                         /* write the data different ways for scalar, vector and tensor */
903                         if (nCompReqd == 1) {
904                           sprintf(tmp_buffer,FLOAT_SCALAR_FORMAT,values[0]);
905                         } else if (nCompReqd == 3) {
906                           if (shape==1) {
907                            sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,values[0],0.,0.);
908                           } else if (shape==2) {
909                            sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,values[0],values[1],0.);
910                           } else if (shape==3) {
911                            sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,values[0],values[1],values[2]);
912                           }
913                         } else if (nCompReqd == 9) {
914                           if (shape==1) {
915                            sprintf(tmp_buffer,FLOAT_TENSOR_FORMAT,values[0],0.,0.,
916                                                                   0.,0.,0.,
917                                                                   0.,0.,0.);
918                           } else if (shape==2) {
919                            sprintf(tmp_buffer,FLOAT_TENSOR_FORMAT,values[0],values[1],0.,
920                                                                   values[2],values[3],0.,
921                                                                   0.,0.,0.);
922                           } else if (shape==3) {
923                            sprintf(tmp_buffer,FLOAT_TENSOR_FORMAT,values[0],values[1],values[2],
924                                                                   values[3],values[4],values[5],
925                                                                   values[6],values[7],values[8]);
926                           }
927                         }
928                         if ( mpi_size > 1) {
929                           __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use);
930                         } else {
931                           fprintf(fileHandle_p,tmp_buffer);
932                         }
933                      }
934                   }
935                   if ( mpi_size > 1) {
936                       #ifdef PASO_MPI
937                         MPI_File_write_ordered(mpi_fileHandle_p,txt_buffer,txt_buffer_in_use, MPI_CHAR, &mpi_status);
938                       #endif    
939                       if ( my_mpi_rank == 0) {
940                          #ifdef PASO_MPI
941                             MPI_File_iwrite_shared(mpi_fileHandle_p,tag_End_DataArray,strlen(tag_End_DataArray),MPI_CHAR,&mpi_req);
942                             MPI_Wait(&mpi_req,&mpi_status);
943                          #endif
944                       }
945                   } else {
946                      fprintf(fileHandle_p,tag_End_DataArray);
947                   }
948                }
949              }
950            }
951            if ( mpi_size > 1) {
952              if ( my_mpi_rank == 0) {
953                 #ifdef PASO_MPI
954                    MPI_File_iwrite_shared(mpi_fileHandle_p,tag_End_PointData,strlen(tag_End_PointData),MPI_CHAR,&mpi_req);
955                    MPI_Wait(&mpi_req,&mpi_status);
956                 #endif
957              }
958            } else {
959                fprintf(fileHandle_p,tag_End_PointData);
960            }
961      }
962      if (Finley_noError()) {
963         if ( mpi_size > 1) {
964           if ( my_mpi_rank == 0) {
965              #ifdef PASO_MPI
966                 MPI_File_iwrite_shared(mpi_fileHandle_p,footer,strlen(footer),MPI_CHAR,&mpi_req);
967                 MPI_Wait(&mpi_req,&mpi_status);
968                 #ifdef MPIO_HINTS
969                   MPI_Info_free(&mpi_info);
970                   #undef MPIO_HINTS
971                 #endif
972              #endif
973            }
974            #ifdef PASO_MPI
975               MPI_File_close(&mpi_fileHandle_p);
976            #endif
977         } else {
978             fprintf(fileHandle_p,footer);
979             fclose(fileHandle_p);
980         }
981      }
982      TMPMEMFREE(isCellCentered);
983      TMPMEMFREE(txt_buffer);
984    return;    return;
985    #else
986      /* Don't kill the job if saveVTK() doesn't work */
987      fprintf(stderr, "\n\nsaveVTK warning: VTK is not available\n\n\n");
988    #endif
989  }  }
   
 /*  
  * $Log$  
  * 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.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.147  
changed lines
  Added in v.1705

  ViewVC Help
Powered by ViewVC 1.1.26