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

Contents of /trunk/dudley/src/Mesh_saveVTK.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3259 - (show annotations)
Mon Oct 11 01:48:14 2010 UTC (9 years ago) by jfenwick
File MIME type: text/plain
File size: 30379 byte(s)
Merging dudley and scons updates from branches

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26