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

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

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

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

Legend:
Removed from v.123  
changed lines
  Added in v.1312

  ViewVC Help
Powered by ViewVC 1.1.26