/[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 3079 - (hide annotations)
Tue Aug 3 04:04:51 2010 UTC (8 years, 9 months ago) by jfenwick
Original Path: branches/domexper/finley/src/Mesh_saveVTK.c
File MIME type: text/plain
File size: 46445 byte(s)
Some experiments on finley

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26