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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2748 - (hide annotations)
Tue Nov 17 07:32:59 2009 UTC (9 years, 9 months ago) by gross
File MIME type: text/plain
File size: 46605 byte(s)
Macro elements are implemented now. VTK writer for macro elements still needs testing.
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 caltinay 2744 #define LEN_INT_FORMAT (unsigned int)(9+2)
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 2744 /* strlen("-1.234567e+789 ") == 15 */
31     #define LEN_TENSOR_FORMAT (unsigned int)(9*15+2)
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 gross 2748 if ( (type == Rec9) || (type == Rec9Macro) ) {
99 caltinay 2141 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 gross 2748 } else if ((type == Hex27) || (type == Hex27Macro) ){
110 caltinay 2141 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 gross 2748 ElementTypeId typeId = NoRef;
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 gross 2748 typeId = elements->referenceElementSet->referenceElement->LinearType->TypeId;
379 caltinay 2141 if (hasReducedElements) {
380 gross 2748 quadNodes_p=elements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->QuadNodes;
381 caltinay 2141 } else {
382 gross 2748 quadNodes_p=elements->referenceElementSet->referenceElement->Parametrization->QuadNodes;
383 caltinay 2141 }
384     } else {
385 gross 2748 typeId = elements->referenceElementSet->referenceElement->Type->TypeId;
386     if (hasReducedElements) {
387     quadNodes_p=elements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->QuadNodes;
388     } else {
389     quadNodes_p=elements->referenceElementSet->referenceElement->Parametrization->QuadNodes;
390     }
391 caltinay 2141 }
392     switch (typeId) {
393     case Point1:
394     case Line2Face:
395     case Line3Face:
396     case Point1_Contact:
397     case Line2Face_Contact:
398     case Line3Face_Contact:
399     cellType = VTK_VERTEX;
400     numVTKNodesPerElement = 1;
401     break;
402 dhawcroft 793
403 caltinay 2141 case Line2:
404     case Tri3Face:
405     case Rec4Face:
406     case Line2_Contact:
407     case Tri3_Contact:
408     case Tri3Face_Contact:
409     case Rec4Face_Contact:
410     cellType = VTK_LINE;
411     numVTKNodesPerElement = 2;
412     break;
413 dhawcroft 793
414 caltinay 2141 case Tri3:
415     case Tet4Face:
416     case Tet4Face_Contact:
417     cellType = VTK_TRIANGLE;
418     numVTKNodesPerElement = 3;
419     break;
420 dhawcroft 793
421 caltinay 2141 case Rec4:
422     case Hex8Face:
423     case Rec4_Contact:
424     case Hex8Face_Contact:
425     cellType = VTK_QUAD;
426     numVTKNodesPerElement = 4;
427     break;
428 dhawcroft 793
429 gross 2748 case Rec9Macro:
430     case Rec9:
431 caltinay 2141 numCellFactor = 4;
432     cellType = VTK_QUAD;
433     numVTKNodesPerElement = 4;
434     break;
435 dhawcroft 793
436 caltinay 2141 case Tet4:
437     cellType = VTK_TETRA;
438     numVTKNodesPerElement = 4;
439     break;
440 dhawcroft 793
441 caltinay 2141 case Hex8:
442     cellType = VTK_HEXAHEDRON;
443     numVTKNodesPerElement = 8;
444     break;
445 dhawcroft 796
446 caltinay 2141 case Line3:
447 gross 2748 case Line3Macro:
448 caltinay 2141 case Tri6Face:
449     case Rec8Face:
450     case Line3_Contact:
451     case Tri6Face_Contact:
452     case Rec8Face_Contact:
453     cellType = VTK_QUADRATIC_EDGE;
454     numVTKNodesPerElement = 3;
455     break;
456 dhawcroft 793
457 caltinay 2141 case Tri6:
458 gross 2748 case Tri6Macro:
459 caltinay 2141 case Tet10Face:
460     case Tri6_Contact:
461     case Tet10Face_Contact:
462     cellType = VTK_QUADRATIC_TRIANGLE;
463     numVTKNodesPerElement = 6;
464     break;
465 dhawcroft 793
466 caltinay 2141 case Rec8:
467     case Hex20Face:
468     case Rec8_Contact:
469     case Hex20Face_Contact:
470     cellType = VTK_QUADRATIC_QUAD;
471     numVTKNodesPerElement = 8;
472     break;
473 dhawcroft 793
474 gross 2748 case Tet10Macro:
475 caltinay 2141 case Tet10:
476     cellType = VTK_QUADRATIC_TETRA;
477     numVTKNodesPerElement = 10;
478     break;
479 gross 1741
480 caltinay 2141 case Hex20:
481     cellType = VTK_QUADRATIC_HEXAHEDRON;
482     numVTKNodesPerElement = 20;
483     break;
484 gross 1741
485 gross 2748 case Hex27Macro:
486 caltinay 2141 case Hex27:
487     numCellFactor = 8;
488     cellType = VTK_HEXAHEDRON;
489     numVTKNodesPerElement = 8;
490     break;
491 gross 1741
492 caltinay 2141 default:
493 gross 2748 sprintf(errorMsg, "saveVTK: Element type %s is not supported by VTK.", elements->referenceElementSet->referenceElement->Type->Name);
494 caltinay 2141 Finley_setError(VALUE_ERROR, errorMsg);
495     }
496     }
497     }
498 gross 1741
499 caltinay 2141 /* allocate enough memory for text buffer */
500 gross 1741
501 gross 2421 txtBufferSize = strlen(vtkHeader) + 3*LEN_INT_FORMAT + (30+3*maxNameLen)+strlen(metadata)+strlen(metadata_schema);
502 caltinay 2141 if (mpi_size > 1) {
503 caltinay 2743 txtBufferSize = MAX(txtBufferSize, myNumPoints * LEN_TMP_BUFFER);
504 caltinay 2141 txtBufferSize = MAX(txtBufferSize, numCellFactor * myNumCells *
505     (LEN_INT_FORMAT * numVTKNodesPerElement + 1));
506     txtBufferSize = MAX(txtBufferSize,
507     numCellFactor * myNumCells * LEN_TENSOR_FORMAT);
508     txtBufferSize = MAX(txtBufferSize, myNumPoints * LEN_TENSOR_FORMAT);
509     }
510     txtBuffer = TMPMEMALLOC(txtBufferSize+1, char);
511 gross 1741
512 caltinay 2141 /* sets error if memory allocation failed */
513     Finley_checkPtr(txtBuffer);
514 gross 1741
515 caltinay 2141 /************************************************************************/
516     /* write number of points and the mesh component */
517 gross 1741
518 caltinay 2141 if (Finley_noError()) {
519     const index_t *nodeIndex;
520     if (FINLEY_REDUCED_NODES == nodeType) {
521 gross 2748 nodeIndex = elements->referenceElementSet->referenceElement->Type->linearNodes;
522     } else if ( (Rec9 == typeId) || (Rec9Macro == typeId) ) {
523 caltinay 2141 nodeIndex = VTK_REC9_INDEX;
524     } else if (Hex20 == typeId) {
525     nodeIndex = VTK_HEX20_INDEX;
526 gross 2748 } else if ( (Hex27 == typeId) || (Hex27Macro == typeId) ){
527 caltinay 2141 nodeIndex = VTK_HEX27_INDEX;
528 gross 2748 } else if (numVTKNodesPerElement != elements->referenceElementSet->referenceElement->Type->numNodes) {
529     nodeIndex = elements->referenceElementSet->referenceElement->Type->relevantGeoNodes;
530 ksteube 1312 } else {
531 caltinay 2141 nodeIndex = NULL;
532 ksteube 1312 }
533 gross 1741
534 gross 2421 if (strlen(metadata)>0) {
535     if (strlen(metadata_schema)>0) {
536     sprintf(txtBuffer, vtkHeader," ",metadata_schema,metadata,"\n",globalNumPoints, numCellFactor*globalNumCells, 3);
537     } else {
538     sprintf(txtBuffer, vtkHeader,"","",metadata,"\n",globalNumPoints, numCellFactor*globalNumCells, 3);
539     }
540     } else {
541     if (strlen(metadata_schema)>0) {
542     sprintf(txtBuffer, vtkHeader," ",metadata_schema,"","",globalNumPoints, numCellFactor*globalNumCells, 3);
543     } else {
544     sprintf(txtBuffer, vtkHeader,"","","","",globalNumPoints, numCellFactor*globalNumCells, 3);
545     }
546     }
547 gross 1741
548 caltinay 2141 if (mpi_size > 1) {
549     /* write the nodes */
550     MPI_RANK0_WRITE_SHARED(txtBuffer);
551     txtBuffer[0] = '\0';
552     txtBufferInUse = 0;
553     if (nDim==2) {
554     for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
555     if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
556     sprintf(tmpBuffer, VECTOR_FORMAT,
557     mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
558     mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
559     0.);
560     __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
561     }
562     }
563     } else {
564     for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
565     if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
566     sprintf(tmpBuffer, VECTOR_FORMAT,
567     mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
568     mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
569     mesh_p->Nodes->Coordinates[INDEX2(2, i, nDim)]);
570     __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
571     }
572     }
573     } /* nDim */
574 gross 2385 MPI_WRITE_ORDERED(txtBuffer);
575 gross 1741
576 caltinay 2141 /* write the cells */
577     MPI_RANK0_WRITE_SHARED(tags_End_Points_and_Start_Conn);
578     txtBuffer[0] = '\0';
579     txtBufferInUse = 0;
580     if (nodeIndex == NULL) {
581     for (i = 0; i < numCells; i++) {
582     if (elements->Owner[i] == my_mpi_rank) {
583     for (j = 0; j < numVTKNodesPerElement; j++) {
584     sprintf(tmpBuffer, INT_FORMAT, globalNodeIndex[elements->Nodes[INDEX2(j, i, NN)]]);
585     __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
586     }
587     __STRCAT(txtBuffer, NEWLINE, txtBufferInUse);
588     }
589     }
590     } else {
591     for (i = 0; i < numCells; i++) {
592     if (elements->Owner[i] == my_mpi_rank) {
593     for (l = 0; l < numCellFactor; l++) {
594     const int* idx=&nodeIndex[l*numVTKNodesPerElement];
595     for (j = 0; j < numVTKNodesPerElement; j++) {
596     sprintf(tmpBuffer, INT_FORMAT, globalNodeIndex[elements->Nodes[INDEX2(idx[j], i, NN)]]);
597     __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
598     }
599     __STRCAT(txtBuffer, NEWLINE, txtBufferInUse);
600     }
601     }
602     }
603     } /* nodeIndex */
604 gross 2385 MPI_WRITE_ORDERED(txtBuffer);
605 gross 1741
606 caltinay 2141 /* write the offsets */
607     MPI_RANK0_WRITE_SHARED(tags_End_Conn_and_Start_Offset);
608     txtBuffer[0] = '\0';
609     txtBufferInUse = 0;
610     for (i = numVTKNodesPerElement*(myFirstCell*numCellFactor+1);
611     i <= (myFirstCell+myNumCells)*numVTKNodesPerElement*numCellFactor; i += numVTKNodesPerElement)
612     {
613     sprintf(tmpBuffer, INT_NEWLINE_FORMAT, i);
614     __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
615     }
616 gross 2385 MPI_WRITE_ORDERED(txtBuffer);
617 gross 1741
618 caltinay 2141 /* write element type */
619     sprintf(tmpBuffer, INT_NEWLINE_FORMAT, cellType);
620     MPI_RANK0_WRITE_SHARED(tags_End_Offset_and_Start_Type);
621     txtBuffer[0] = '\0';
622     txtBufferInUse = 0;
623     for (i = numVTKNodesPerElement*(myFirstCell*numCellFactor+1);
624     i <= (myFirstCell+myNumCells)*numVTKNodesPerElement*numCellFactor; i += numVTKNodesPerElement)
625     {
626     __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
627     }
628 gross 2385 MPI_WRITE_ORDERED(txtBuffer);
629 caltinay 2141 /* finalize cell information */
630     strcpy(txtBuffer, "</DataArray>\n</Cells>\n");
631     MPI_RANK0_WRITE_SHARED(txtBuffer);
632 gross 1741
633 caltinay 2141 } else { /***** mpi_size == 1 *****/
634 gross 1741
635 caltinay 2141 /* write the nodes */
636     fputs(txtBuffer, fileHandle_p);
637     if (nDim==2) {
638     for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
639     if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
640     fprintf(fileHandle_p, VECTOR_FORMAT,
641     mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
642     mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
643     0.);
644     }
645     }
646     } else {
647     for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
648     if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
649     fprintf(fileHandle_p, VECTOR_FORMAT,
650     mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
651     mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
652     mesh_p->Nodes->Coordinates[INDEX2(2, i, nDim)]);
653     }
654     }
655     } /* nDim */
656 gross 1741
657 caltinay 2141 /* write the cells */
658     fputs(tags_End_Points_and_Start_Conn, fileHandle_p);
659     if (nodeIndex == NULL) {
660     for (i = 0; i < numCells; i++) {
661     for (j = 0; j < numVTKNodesPerElement; j++) {
662     fprintf(fileHandle_p, INT_FORMAT, globalNodeIndex[elements->Nodes[INDEX2(j, i, NN)]]);
663     }
664     fprintf(fileHandle_p, NEWLINE);
665     }
666     } else {
667     for (i = 0; i < numCells; i++) {
668     for (l = 0; l < numCellFactor; l++) {
669     const int* idx=&nodeIndex[l*numVTKNodesPerElement];
670     for (j = 0; j < numVTKNodesPerElement; j++) {
671     fprintf(fileHandle_p, INT_FORMAT, globalNodeIndex[elements->Nodes[INDEX2(idx[j], i, NN)]]);
672     }
673     fprintf(fileHandle_p, NEWLINE);
674     }
675     }
676     } /* nodeIndex */
677 gross 1741
678 caltinay 2141 /* write the offsets */
679     fputs(tags_End_Conn_and_Start_Offset, fileHandle_p);
680     for (i = numVTKNodesPerElement; i <= numCells*numVTKNodesPerElement*numCellFactor; i += numVTKNodesPerElement) {
681     fprintf(fileHandle_p, INT_NEWLINE_FORMAT, i);
682     }
683 gross 1741
684 caltinay 2141 /* write element type */
685     sprintf(tmpBuffer, INT_NEWLINE_FORMAT, cellType);
686     fputs(tags_End_Offset_and_Start_Type, fileHandle_p);
687     for (i = 0; i < numCells*numCellFactor; i++)
688     fputs(tmpBuffer, fileHandle_p);
689     /* finalize cell information */
690     fputs("</DataArray>\n</Cells>\n", fileHandle_p);
691     } /* mpi_size */
692 gross 1741
693 caltinay 2141 } /* Finley_noError */
694 dhawcroft 796
695 caltinay 2141 /************************************************************************/
696     /* write cell data */
697 dhawcroft 796
698 caltinay 2141 if (writeCellData && Finley_noError()) {
699     bool_t set_scalar=FALSE, set_vector=FALSE, set_tensor=FALSE;
700     /* mark the active data arrays */
701     strcpy(txtBuffer, "<CellData");
702     for (dataIdx = 0; dataIdx < num_data; dataIdx++) {
703     if (!isEmpty(data_pp[dataIdx]) && isCellCentered[dataIdx]) {
704     /* rank == 0 <--> scalar data */
705     /* rank == 1 <--> vector data */
706     /* rank == 2 <--> tensor data */
707     switch (getDataPointRank(data_pp[dataIdx])) {
708     case 0:
709     if (!set_scalar) {
710     strcat(txtBuffer, " Scalars=\"");
711     strcat(txtBuffer, names_p[dataIdx]);
712     strcat(txtBuffer, "\"");
713     set_scalar = TRUE;
714     }
715     break;
716     case 1:
717     if (!set_vector) {
718     strcat(txtBuffer, " Vectors=\"");
719     strcat(txtBuffer, names_p[dataIdx]);
720     strcat(txtBuffer, "\"");
721     set_vector = TRUE;
722     }
723     break;
724     case 2:
725     if (!set_tensor) {
726     strcat(txtBuffer, " Tensors=\"");
727     strcat(txtBuffer, names_p[dataIdx]);
728     strcat(txtBuffer, "\"");
729     set_tensor = TRUE;
730     }
731     break;
732     default:
733     sprintf(errorMsg, "saveVTK: data %s: VTK supports data with rank <= 2 only.", names_p[dataIdx]);
734     Finley_setError(VALUE_ERROR, errorMsg);
735     }
736     }
737     if (!Finley_noError())
738     break;
739 dhawcroft 818 }
740 dhawcroft 793 }
741 caltinay 2141 /* only continue if no error occurred */
742     if (writeCellData && Finley_noError()) {
743     strcat(txtBuffer, ">\n");
744     if ( mpi_size > 1) {
745     MPI_RANK0_WRITE_SHARED(txtBuffer);
746     } else {
747     fputs(txtBuffer, fileHandle_p);
748 dhawcroft 793 }
749 dhawcroft 818
750 caltinay 2141 /* write the arrays */
751     for (dataIdx = 0; dataIdx < num_data; dataIdx++) {
752     if (!isEmpty(data_pp[dataIdx]) && isCellCentered[dataIdx]) {
753     dim_t numPointsPerSample=getNumDataPointsPerSample(data_pp[dataIdx]);
754     dim_t rank = getDataPointRank(data_pp[dataIdx]);
755     dim_t nComp = getDataPointSize(data_pp[dataIdx]);
756     dim_t nCompReqd = 1; /* number of components mpi_required */
757     if (rank == 0) {
758     nCompReqd = 1;
759     shape = 0;
760     } else if (rank == 1) {
761     nCompReqd = 3;
762     shape = getDataPointShape(data_pp[dataIdx], 0);
763     if (shape > 3) {
764     Finley_setError(VALUE_ERROR, "saveVTK: rank 1 objects must have 3 components at most.");
765     }
766     } else {
767     nCompReqd = 9;
768     shape = getDataPointShape(data_pp[dataIdx], 0);
769     if (shape > 3 || shape != getDataPointShape(data_pp[dataIdx], 1)) {
770     Finley_setError(VALUE_ERROR, "saveVTK: rank 2 objects of shape 2x2 or 3x3 supported only.");
771     }
772     }
773     /* bail out if an error occurred */
774     if (!Finley_noError())
775     break;
776 dhawcroft 793
777 caltinay 2141 sprintf(txtBuffer, tag_Float_DataArray, names_p[dataIdx], nCompReqd);
778     if ( mpi_size > 1) {
779     MPI_RANK0_WRITE_SHARED(txtBuffer);
780     } else {
781     fputs(txtBuffer, fileHandle_p);
782     }
783 gross 1741
784 caltinay 2141 txtBuffer[0] = '\0';
785     txtBufferInUse = 0;
786     for (i=0; i<numCells; i++) {
787     if (elements->Owner[i] == my_mpi_rank) {
788 caltinay 2743 void* sampleBuffer=allocSampleBuffer(data_pp[dataIdx]);
789 jfenwick 2271 __const double *values = getSampleDataRO(data_pp[dataIdx], i,sampleBuffer);
790 caltinay 2141 for (l = 0; l < numCellFactor; l++) {
791     double sampleAvg[NCOMP_MAX];
792     dim_t nCompUsed = MIN(nComp, NCOMP_MAX);
793    
794     /* average over number of points in the sample */
795     if (isExpanded(data_pp[dataIdx])) {
796 gross 2748 dim_t hits=0;
797 caltinay 2141 for (k=0; k<nCompUsed; k++) sampleAvg[k]=0;
798     for (j=0; j<numPointsPerSample; j++) {
799     if (nodeInQuadrant(quadNodes_p, typeId, j, l)) {
800     hits++;
801     for (k=0; k<nCompUsed; k++) {
802     sampleAvg[k] += values[INDEX2(k,j,nComp)];
803     }
804 gross 1741 }
805 caltinay 2141 }
806     for (k=0; k<nCompUsed; k++)
807     sampleAvg[k] /= MAX(hits, 1);
808     } else {
809     for (k=0; k<nCompUsed; k++)
810     sampleAvg[k] = values[k];
811     } /* isExpanded */
812    
813     /* if the number of required components is more than
814     * the number of actual components, pad with zeros
815     */
816     /* probably only need to get shape of first element */
817     if (nCompReqd == 1) {
818     sprintf(tmpBuffer, SCALAR_FORMAT, sampleAvg[0]);
819     } else if (nCompReqd == 3) {
820     if (shape==1) {
821     sprintf(tmpBuffer, VECTOR_FORMAT,
822     sampleAvg[0], 0.f, 0.f);
823     } else if (shape==2) {
824     sprintf(tmpBuffer, VECTOR_FORMAT,
825     sampleAvg[0], sampleAvg[1], 0.f);
826     } else if (shape==3) {
827     sprintf(tmpBuffer, VECTOR_FORMAT,
828     sampleAvg[0],sampleAvg[1],sampleAvg[2]);
829     }
830     } else if (nCompReqd == 9) {
831     if (shape==1) {
832     sprintf(tmpBuffer, TENSOR_FORMAT,
833     sampleAvg[0], 0.f, 0.f,
834     0.f, 0.f, 0.f,
835     0.f, 0.f, 0.f);
836     } else if (shape==2) {
837     sprintf(tmpBuffer, TENSOR_FORMAT,
838     sampleAvg[0], sampleAvg[1], 0.f,
839     sampleAvg[2], sampleAvg[3], 0.f,
840     0.f, 0.f, 0.f);
841     } else if (shape==3) {
842     sprintf(tmpBuffer, TENSOR_FORMAT,
843     sampleAvg[0],sampleAvg[1],sampleAvg[2],
844     sampleAvg[3],sampleAvg[4],sampleAvg[5],
845     sampleAvg[6],sampleAvg[7],sampleAvg[8]);
846     }
847     }
848 gross 1741 if ( mpi_size > 1) {
849 caltinay 2141 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
850 gross 1741 } else {
851 caltinay 2141 fputs(tmpBuffer, fileHandle_p);
852 gross 1741 }
853 caltinay 2141 } /* for l (numCellFactor) */
854 caltinay 2743 freeSampleBuffer(sampleBuffer);
855 caltinay 2141 } /* if I am the owner */
856     } /* for i (numCells) */
857    
858     if ( mpi_size > 1) {
859 gross 2385 MPI_WRITE_ORDERED(txtBuffer);
860 caltinay 2141 MPI_RANK0_WRITE_SHARED(tag_End_DataArray);
861     } else {
862     fputs(tag_End_DataArray, fileHandle_p);
863     }
864     } /* !isEmpty && cellCentered */
865     } /* for dataIdx */
866    
867     strcpy(txtBuffer, "</CellData>\n");
868     if ( mpi_size > 1) {
869     MPI_RANK0_WRITE_SHARED(txtBuffer);
870     } else {
871     fputs(txtBuffer, fileHandle_p);
872 dhawcroft 793 }
873 caltinay 2141 } /* if noError && writeCellData */
874 dhawcroft 793
875 caltinay 2141 /************************************************************************/
876     /* write point data */
877    
878     if (writePointData && Finley_noError()) {
879     /* mark the active data arrays */
880     bool_t set_scalar=FALSE, set_vector=FALSE, set_tensor=FALSE;
881     strcpy(txtBuffer, "<PointData");
882     for (dataIdx = 0; dataIdx<num_data; dataIdx++) {
883     if (!isEmpty(data_pp[dataIdx]) && !isCellCentered[dataIdx]) {
884     switch (getDataPointRank(data_pp[dataIdx])) {
885     case 0:
886     if (!set_scalar) {
887     strcat(txtBuffer, " Scalars=\"");
888     strcat(txtBuffer, names_p[dataIdx]);
889     strcat(txtBuffer, "\"");
890     set_scalar = TRUE;
891     }
892     break;
893     case 1:
894     if (!set_vector) {
895     strcat(txtBuffer, " Vectors=\"");
896     strcat(txtBuffer, names_p[dataIdx]);
897     strcat(txtBuffer, "\"");
898     set_vector = TRUE;
899     }
900     break;
901     case 2:
902     if (!set_tensor) {
903     strcat(txtBuffer, " Tensors=\"");
904     strcat(txtBuffer, names_p[dataIdx]);
905     strcat(txtBuffer, "\"");
906     set_tensor = TRUE;
907     }
908     break;
909     default:
910     sprintf(errorMsg, "saveVTK: data %s: VTK supports data with rank <= 2 only.", names_p[dataIdx]);
911     Finley_setError(VALUE_ERROR, errorMsg);
912     }
913 dhawcroft 793 }
914 caltinay 2141 if (!Finley_noError())
915     break;
916 dhawcroft 793 }
917 caltinay 2141 }
918     /* only continue if no error occurred */
919     if (writePointData && Finley_noError()) {
920     strcat(txtBuffer, ">\n");
921     if ( mpi_size > 1) {
922     MPI_RANK0_WRITE_SHARED(txtBuffer);
923     } else {
924     fputs(txtBuffer, fileHandle_p);
925 gross 1062 }
926 caltinay 2141
927     /* write the arrays */
928     for (dataIdx = 0; dataIdx < num_data; dataIdx++) {
929     if (!isEmpty(data_pp[dataIdx]) && !isCellCentered[dataIdx]) {
930     Finley_NodeMapping* nodeMapping;
931     dim_t rank = getDataPointRank(data_pp[dataIdx]);
932     dim_t nCompReqd = 1; /* number of components mpi_required */
933     if (getFunctionSpaceType(data_pp[dataIdx]) == FINLEY_REDUCED_NODES) {
934     nodeMapping = mesh_p->Nodes->reducedNodesMapping;
935     } else {
936     nodeMapping = mesh_p->Nodes->nodesMapping;
937     }
938     if (rank == 0) {
939     nCompReqd = 1;
940     shape = 0;
941     } else if (rank == 1) {
942     nCompReqd = 3;
943     shape = getDataPointShape(data_pp[dataIdx], 0);
944     if (shape > 3) {
945     Finley_setError(VALUE_ERROR, "saveVTK: rank 1 objects must have 3 components at most.");
946     }
947     } else {
948     nCompReqd = 9;
949     shape=getDataPointShape(data_pp[dataIdx], 0);
950     if (shape > 3 || shape != getDataPointShape(data_pp[dataIdx], 1)) {
951     Finley_setError(VALUE_ERROR, "saveVTK: rank 2 objects of shape 2x2 or 3x3 supported only.");
952     }
953     }
954     /* bail out if an error occurred */
955     if (!Finley_noError())
956     break;
957    
958     sprintf(txtBuffer, tag_Float_DataArray, names_p[dataIdx], nCompReqd);
959     if ( mpi_size > 1) {
960     MPI_RANK0_WRITE_SHARED(txtBuffer);
961     } else {
962     fputs(txtBuffer, fileHandle_p);
963     }
964    
965     txtBuffer[0] = '\0';
966     txtBufferInUse = 0;
967     for (i=0; i<mesh_p->Nodes->numNodes; i++) {
968     k = globalNodeIndex[i];
969     if ( (myFirstNode <= k) && (k < myLastNode) ) {
970 caltinay 2743 void* sampleBuffer=allocSampleBuffer(data_pp[dataIdx]);
971 jfenwick 2271 __const double *values = getSampleDataRO(data_pp[dataIdx], nodeMapping->target[i], sampleBuffer);
972 caltinay 2141 /* if the number of mpi_required components is more than
973     * the number of actual components, pad with zeros.
974     * Probably only need to get shape of first element */
975     if (nCompReqd == 1) {
976     sprintf(tmpBuffer, SCALAR_FORMAT, values[0]);
977     } else if (nCompReqd == 3) {
978     if (shape==1) {
979     sprintf(tmpBuffer, VECTOR_FORMAT,
980     values[0], 0.f, 0.f);
981     } else if (shape==2) {
982     sprintf(tmpBuffer, VECTOR_FORMAT,
983     values[0], values[1], 0.f);
984     } else if (shape==3) {
985     sprintf(tmpBuffer, VECTOR_FORMAT,
986     values[0], values[1], values[2]);
987     }
988     } else if (nCompReqd == 9) {
989     if (shape==1) {
990     sprintf(tmpBuffer, TENSOR_FORMAT,
991     values[0], 0.f, 0.f,
992     0.f, 0.f, 0.f,
993     0.f, 0.f, 0.f);
994     } else if (shape==2) {
995     sprintf(tmpBuffer, TENSOR_FORMAT,
996     values[0], values[1], 0.f,
997     values[2], values[3], 0.f,
998     0.f, 0.f, 0.f);
999     } else if (shape==3) {
1000     sprintf(tmpBuffer, TENSOR_FORMAT,
1001     values[0], values[1], values[2],
1002     values[3], values[4], values[5],
1003     values[6], values[7], values[8]);
1004     }
1005     }
1006     if ( mpi_size > 1) {
1007     __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
1008     } else {
1009     fputs(tmpBuffer, fileHandle_p);
1010     }
1011 caltinay 2743 freeSampleBuffer(sampleBuffer); /* no-one needs values anymore */
1012 caltinay 2141 } /* if this is my node */
1013     } /* for i (numNodes) */
1014    
1015     if ( mpi_size > 1) {
1016 gross 2385 MPI_WRITE_ORDERED(txtBuffer);
1017 caltinay 2141 MPI_RANK0_WRITE_SHARED(tag_End_DataArray);
1018     } else {
1019     fputs(tag_End_DataArray, fileHandle_p);
1020     }
1021     } /* !isEmpty && !isCellCentered */
1022     } /* for dataIdx */
1023    
1024     strcpy(txtBuffer, "</PointData>\n");
1025 ksteube 1312 if ( mpi_size > 1) {
1026 caltinay 2141 MPI_RANK0_WRITE_SHARED(txtBuffer);
1027 gross 1062 } else {
1028 caltinay 2141 fputs(txtBuffer, fileHandle_p);
1029 gross 1062 }
1030 caltinay 2141 } /* if noError && writePointData */
1031    
1032     /* Final write to VTK file */
1033     if (Finley_noError()) {
1034     if (mpi_size > 1) {
1035     MPI_RANK0_WRITE_SHARED(vtkFooter);
1036     } else {
1037     fputs(vtkFooter, fileHandle_p);
1038 dhawcroft 793 }
1039 caltinay 2141 }
1040    
1041     if ( mpi_size > 1) {
1042 caltinay 2094 #ifdef PASO_MPI
1043 caltinay 2141 MPI_File_close(&mpi_fileHandle_p);
1044 caltinay 2743 MPI_Barrier(mesh_p->Nodes->MPIInfo->comm);
1045 caltinay 2094 #endif
1046 caltinay 2141 } else {
1047     fclose(fileHandle_p);
1048     }
1049     TMPMEMFREE(isCellCentered);
1050     TMPMEMFREE(txtBuffer);
1051 jgs 110 }
1052 caltinay 2141

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26