/[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 3227 - (hide annotations)
Thu Sep 30 06:07:08 2010 UTC (8 years, 8 months ago) by jfenwick
File MIME type: text/plain
File size: 30918 byte(s)
Pass1 or moving MPI stuff out of paso

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26