/[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 3086 - (hide annotations)
Thu Aug 5 05:07:58 2010 UTC (8 years, 9 months ago) by jfenwick
File MIME type: text/plain
File size: 46445 byte(s)
Another pass at removing 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 jfenwick 3086 /* Nodal data needs to be given on DUDLEY_NODES or DUDLEY_REDUCED_NODES */
17 caltinay 2141 /***************************************************************************/
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 jfenwick 3086 void Dudley_Mesh_saveVTK(const char *filename_p,
150     Dudley_Mesh *mesh_p,
151 ksteube 1312 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 jfenwick 3086 int cellType=0, nodeType=DUDLEY_NODES, elementType=DUDLEY_UNKNOWN;
179     Dudley_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 jfenwick 3086 Dudley_setError(TYPE_ERROR, "saveVTK: spatial dimension 2 or 3 is supported only.");
220 caltinay 2141 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 jfenwick 3086 Dudley_setError(IO_ERROR, errorMsg);
240 caltinay 2141 } 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 jfenwick 3086 Dudley_setError(IO_ERROR, errorMsg);
250 caltinay 2141 }
251     }
252     if (!Paso_MPIInfo_noError(mesh_p->Nodes->MPIInfo)) return;
253 dhawcroft 793
254 jfenwick 3086 /* General note: From this point if an error occurs Dudley_setError is
255 caltinay 2141 * 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 jfenwick 3086 if (!Dudley_checkPtr(isCellCentered)) {
264 caltinay 2141 for (dataIdx=0; dataIdx<num_data; ++dataIdx) {
265     if (! isEmpty(data_pp[dataIdx])) {
266     switch(getFunctionSpaceType(data_pp[dataIdx]) ) {
267 jfenwick 3086 case DUDLEY_NODES:
268     nodeType = (nodeType == DUDLEY_REDUCED_NODES) ? DUDLEY_REDUCED_NODES : DUDLEY_NODES;
269 caltinay 2141 isCellCentered[dataIdx] = FALSE;
270 jfenwick 3086 if (elementType==DUDLEY_UNKNOWN || elementType==DUDLEY_ELEMENTS) {
271     elementType = DUDLEY_ELEMENTS;
272 caltinay 2141 } else {
273 jfenwick 3086 Dudley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
274 caltinay 2141 }
275     break;
276 jfenwick 3086 case DUDLEY_REDUCED_NODES:
277     nodeType = DUDLEY_REDUCED_NODES;
278 caltinay 2141 isCellCentered[dataIdx] = FALSE;
279 jfenwick 3086 if (elementType==DUDLEY_UNKNOWN || elementType==DUDLEY_ELEMENTS) {
280     elementType = DUDLEY_ELEMENTS;
281 caltinay 2141 } else {
282 jfenwick 3086 Dudley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
283 caltinay 2141 }
284     break;
285 jfenwick 3086 case DUDLEY_REDUCED_ELEMENTS:
286 caltinay 2141 hasReducedElements = TRUE;
287 jfenwick 3086 case DUDLEY_ELEMENTS:
288 caltinay 2141 isCellCentered[dataIdx] = TRUE;
289 jfenwick 3086 if (elementType==DUDLEY_UNKNOWN || elementType==DUDLEY_ELEMENTS) {
290     elementType = DUDLEY_ELEMENTS;
291 caltinay 2141 } else {
292 jfenwick 3086 Dudley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
293 caltinay 2141 }
294     break;
295 jfenwick 3086 case DUDLEY_REDUCED_FACE_ELEMENTS:
296 caltinay 2141 hasReducedElements = TRUE;
297 jfenwick 3086 case DUDLEY_FACE_ELEMENTS:
298 caltinay 2141 isCellCentered[dataIdx] = TRUE;
299 jfenwick 3086 if (elementType==DUDLEY_UNKNOWN || elementType==DUDLEY_FACE_ELEMENTS) {
300     elementType = DUDLEY_FACE_ELEMENTS;
301 caltinay 2141 } else {
302 jfenwick 3086 Dudley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
303 caltinay 2141 }
304     break;
305 jfenwick 3086 case DUDLEY_POINTS:
306 caltinay 2141 isCellCentered[dataIdx]=TRUE;
307 jfenwick 3086 if (elementType==DUDLEY_UNKNOWN || elementType==DUDLEY_POINTS) {
308     elementType = DUDLEY_POINTS;
309 caltinay 2141 } else {
310 jfenwick 3086 Dudley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
311 caltinay 2141 }
312     break;
313 jfenwick 3086 case DUDLEY_REDUCED_CONTACT_ELEMENTS_1:
314 caltinay 2141 hasReducedElements = TRUE;
315 jfenwick 3086 case DUDLEY_CONTACT_ELEMENTS_1:
316 caltinay 2141 isCellCentered[dataIdx] = TRUE;
317 jfenwick 3086 if (elementType==DUDLEY_UNKNOWN || elementType==DUDLEY_CONTACT_ELEMENTS_1) {
318     elementType = DUDLEY_CONTACT_ELEMENTS_1;
319 caltinay 2141 } else {
320 jfenwick 3086 Dudley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
321 caltinay 2141 }
322     break;
323 jfenwick 3086 case DUDLEY_REDUCED_CONTACT_ELEMENTS_2:
324 caltinay 2141 hasReducedElements = TRUE;
325 jfenwick 3086 case DUDLEY_CONTACT_ELEMENTS_2:
326 caltinay 2141 isCellCentered[dataIdx] = TRUE;
327 jfenwick 3086 if (elementType==DUDLEY_UNKNOWN || elementType==DUDLEY_CONTACT_ELEMENTS_1) {
328     elementType = DUDLEY_CONTACT_ELEMENTS_1;
329 caltinay 2141 } else {
330 jfenwick 3086 Dudley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
331 caltinay 2141 }
332     break;
333     default:
334     sprintf(errorMsg, "saveVTK: unknown function space type %d",getFunctionSpaceType(data_pp[dataIdx]));
335 jfenwick 3086 Dudley_setError(TYPE_ERROR, errorMsg);
336 caltinay 2141 }
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 jfenwick 3086 if (Dudley_noError()) {
351     if (nodeType == DUDLEY_REDUCED_NODES) {
352     myFirstNode = Dudley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
353     myLastNode = Dudley_NodeFile_getLastReducedNode(mesh_p->Nodes);
354     globalNumPoints = Dudley_NodeFile_getGlobalNumReducedNodes(mesh_p->Nodes);
355     globalNodeIndex = Dudley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
356 caltinay 2141 } else {
357 jfenwick 3086 myFirstNode = Dudley_NodeFile_getFirstNode(mesh_p->Nodes);
358     myLastNode = Dudley_NodeFile_getLastNode(mesh_p->Nodes);
359     globalNumPoints = Dudley_NodeFile_getGlobalNumNodes(mesh_p->Nodes);
360     globalNodeIndex = Dudley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
361 ksteube 1312 }
362 caltinay 2141 myNumPoints = myLastNode - myFirstNode;
363 jfenwick 3086 if (elementType==DUDLEY_UNKNOWN) elementType=DUDLEY_ELEMENTS;
364 caltinay 2141 switch(elementType) {
365 jfenwick 3086 case DUDLEY_ELEMENTS:
366 caltinay 2141 elements = mesh_p->Elements;
367     break;
368 jfenwick 3086 case DUDLEY_FACE_ELEMENTS:
369 caltinay 2141 elements = mesh_p->FaceElements;
370     break;
371 jfenwick 3086 case DUDLEY_POINTS:
372 caltinay 2141 elements = mesh_p->Points;
373     break;
374 jfenwick 3086 case DUDLEY_CONTACT_ELEMENTS_1:
375 caltinay 2141 elements = mesh_p->ContactElements;
376     break;
377     }
378     if (elements==NULL) {
379 jfenwick 3086 Dudley_setError(SYSTEM_ERROR, "saveVTK: undefined element file");
380 caltinay 2141 } else {
381 jfenwick 3086 /* map dudley element type to VTK element type */
382 caltinay 2141 numCells = elements->numElements;
383 jfenwick 3086 globalNumCells = Dudley_ElementFile_getGlobalNumElements(elements);
384     myNumCells = Dudley_ElementFile_getMyNumElements(elements);
385     myFirstCell = Dudley_ElementFile_getFirstElement(elements);
386 caltinay 2141 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 jfenwick 3086 if (nodeType==DUDLEY_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 jfenwick 3086 Dudley_setError(VALUE_ERROR, errorMsg);
506 caltinay 2141 }
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 jfenwick 3086 Dudley_checkPtr(txtBuffer);
525 gross 1741
526 caltinay 2141 /************************************************************************/
527     /* write number of points and the mesh component */
528 gross 1741
529 jfenwick 3086 if (Dudley_noError()) {
530 caltinay 2141 const index_t *nodeIndex;
531 jfenwick 3086 if (DUDLEY_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 jfenwick 3086 } /* Dudley_noError */
707 dhawcroft 796
708 caltinay 2141 /************************************************************************/
709     /* write cell data */
710 dhawcroft 796
711 jfenwick 3086 if (writeCellData && Dudley_noError()) {
712 caltinay 2141 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 jfenwick 3086 Dudley_setError(VALUE_ERROR, errorMsg);
748 caltinay 2141 }
749     }
750 jfenwick 3086 if (!Dudley_noError())
751 caltinay 2141 break;
752 dhawcroft 818 }
753 dhawcroft 793 }
754 caltinay 2141 /* only continue if no error occurred */
755 jfenwick 3086 if (writeCellData && Dudley_noError()) {
756 caltinay 2141 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 jfenwick 3086 Dudley_setError(VALUE_ERROR, "saveVTK: rank 1 objects must have 3 components at most.");
778 caltinay 2141 }
779     } else {
780     nCompReqd = 9;
781     shape = getDataPointShape(data_pp[dataIdx], 0);
782     if (shape > 3 || shape != getDataPointShape(data_pp[dataIdx], 1)) {
783 jfenwick 3086 Dudley_setError(VALUE_ERROR, "saveVTK: rank 2 objects of shape 2x2 or 3x3 supported only.");
784 caltinay 2141 }
785     }
786     /* bail out if an error occurred */
787 jfenwick 3086 if (!Dudley_noError())
788 caltinay 2141 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 jfenwick 3086 if (writePointData && Dudley_noError()) {
890 caltinay 2141 /* 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 jfenwick 3086 Dudley_setError(VALUE_ERROR, errorMsg);
923 caltinay 2141 }
924 dhawcroft 793 }
925 jfenwick 3086 if (!Dudley_noError())
926 caltinay 2141 break;
927 dhawcroft 793 }
928 caltinay 2141 }
929     /* only continue if no error occurred */
930 jfenwick 3086 if (writePointData && Dudley_noError()) {
931 caltinay 2141 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 jfenwick 3086 Dudley_NodeMapping* nodeMapping;
942 caltinay 2141 dim_t rank = getDataPointRank(data_pp[dataIdx]);
943     dim_t nCompReqd = 1; /* number of components mpi_required */
944 jfenwick 3086 if (getFunctionSpaceType(data_pp[dataIdx]) == DUDLEY_REDUCED_NODES) {
945 caltinay 2141 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 jfenwick 3086 Dudley_setError(VALUE_ERROR, "saveVTK: rank 1 objects must have 3 components at most.");
957 caltinay 2141 }
958     } else {
959     nCompReqd = 9;
960     shape=getDataPointShape(data_pp[dataIdx], 0);
961     if (shape > 3 || shape != getDataPointShape(data_pp[dataIdx], 1)) {
962 jfenwick 3086 Dudley_setError(VALUE_ERROR, "saveVTK: rank 2 objects of shape 2x2 or 3x3 supported only.");
963 caltinay 2141 }
964     }
965     /* bail out if an error occurred */
966 jfenwick 3086 if (!Dudley_noError())
967 caltinay 2141 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 jfenwick 3086 if (Dudley_noError()) {
1043 caltinay 2141 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