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

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

Legend:
Removed from v.1384  
changed lines
  Added in v.2754

  ViewVC Help
Powered by ViewVC 1.1.26