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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2141 - (show annotations)
Tue Dec 9 04:32:32 2008 UTC (10 years, 8 months ago) by caltinay
File MIME type: text/plain
File size: 44850 byte(s)
saveVTK: Reviewed and changed parts of this function due to several issues.
Confirmed that it's still working using finley/test/python/generate_dumps.py
(manually inserted saveVTK) with and without MPI on 1 and 4 processors.

Changes include:
- removed unnecessary USE_VTK check - we are not using any VTK functions
- fixed writing of Rec9 & Hex27 connectivity data in MPI mode
- fixed some cases where the file was not closed and memory was not freed
- replaced many instances of fprintf by fputs (simple strings)
- some minor optimizations by moving or initializing variables and combining
  if-clauses
- got rid of ~300 lines of code and switched to consistent indentation making
  the file a bit better readable


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26