/[escript]/branches/domexper/dudley/src/Mesh_saveVTK.c
ViewVC logotype

Annotation of /branches/domexper/dudley/src/Mesh_saveVTK.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2743 - (hide annotations)
Mon Nov 16 05:28:47 2009 UTC (9 years, 6 months ago) by caltinay
Original Path: trunk/finley/src/Mesh_saveVTK.c
File MIME type: text/plain
File size: 46122 byte(s)
Corrected memory allocation size for tensors in saveVTK. Fixes issue 460.

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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26