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

Legend:
Removed from v.113  
changed lines
  Added in v.2141

  ViewVC Help
Powered by ViewVC 1.1.26