/[escript]/branches/domexper/dudley/src/Mesh_saveVTK.c
ViewVC logotype

Contents of /branches/domexper/dudley/src/Mesh_saveVTK.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2271 - (show annotations)
Mon Feb 16 05:08:29 2009 UTC (10 years, 3 months ago) by jfenwick
Original Path: trunk/finley/src/Mesh_saveVTK.c
File MIME type: text/plain
File size: 45206 byte(s)
Merging version 2269 to trunk

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 (unsigned int)(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 (unsigned int)(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 (unsigned int)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 void* sampleBuffer=allocSampleBuffer(data_pp[dataIdx]);
766 __const double *values = getSampleDataRO(data_pp[dataIdx], i,sampleBuffer);
767 for (l = 0; l < numCellFactor; l++) {
768 double sampleAvg[NCOMP_MAX];
769 dim_t nCompUsed = MIN(nComp, NCOMP_MAX);
770
771 /* average over number of points in the sample */
772 if (isExpanded(data_pp[dataIdx])) {
773 dim_t hits=0, hits_old;
774 for (k=0; k<nCompUsed; k++) sampleAvg[k]=0;
775 for (j=0; j<numPointsPerSample; j++) {
776 hits_old=hits;
777 if (nodeInQuadrant(quadNodes_p, typeId, j, l)) {
778 hits++;
779 for (k=0; k<nCompUsed; k++) {
780 sampleAvg[k] += values[INDEX2(k,j,nComp)];
781 }
782 }
783 }
784 for (k=0; k<nCompUsed; k++)
785 sampleAvg[k] /= MAX(hits, 1);
786 } else {
787 for (k=0; k<nCompUsed; k++)
788 sampleAvg[k] = values[k];
789 } /* isExpanded */
790
791 /* if the number of required components is more than
792 * the number of actual components, pad with zeros
793 */
794 /* probably only need to get shape of first element */
795 if (nCompReqd == 1) {
796 sprintf(tmpBuffer, SCALAR_FORMAT, sampleAvg[0]);
797 } else if (nCompReqd == 3) {
798 if (shape==1) {
799 sprintf(tmpBuffer, VECTOR_FORMAT,
800 sampleAvg[0], 0.f, 0.f);
801 } else if (shape==2) {
802 sprintf(tmpBuffer, VECTOR_FORMAT,
803 sampleAvg[0], sampleAvg[1], 0.f);
804 } else if (shape==3) {
805 sprintf(tmpBuffer, VECTOR_FORMAT,
806 sampleAvg[0],sampleAvg[1],sampleAvg[2]);
807 }
808 } else if (nCompReqd == 9) {
809 if (shape==1) {
810 sprintf(tmpBuffer, TENSOR_FORMAT,
811 sampleAvg[0], 0.f, 0.f,
812 0.f, 0.f, 0.f,
813 0.f, 0.f, 0.f);
814 } else if (shape==2) {
815 sprintf(tmpBuffer, TENSOR_FORMAT,
816 sampleAvg[0], sampleAvg[1], 0.f,
817 sampleAvg[2], sampleAvg[3], 0.f,
818 0.f, 0.f, 0.f);
819 } else if (shape==3) {
820 sprintf(tmpBuffer, TENSOR_FORMAT,
821 sampleAvg[0],sampleAvg[1],sampleAvg[2],
822 sampleAvg[3],sampleAvg[4],sampleAvg[5],
823 sampleAvg[6],sampleAvg[7],sampleAvg[8]);
824 }
825 }
826 if ( mpi_size > 1) {
827 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
828 } else {
829 fputs(tmpBuffer, fileHandle_p);
830 }
831 } /* for l (numCellFactor) */
832 freeSampleBuffer(sampleBuffer);
833 } /* if I am the owner */
834 } /* for i (numCells) */
835
836 if ( mpi_size > 1) {
837 MPI_WRITE_ORDERED(txtBuffer, txtBufferInUse);
838 MPI_RANK0_WRITE_SHARED(tag_End_DataArray);
839 } else {
840 fputs(tag_End_DataArray, fileHandle_p);
841 }
842 } /* !isEmpty && cellCentered */
843 } /* for dataIdx */
844
845 strcpy(txtBuffer, "</CellData>\n");
846 if ( mpi_size > 1) {
847 MPI_RANK0_WRITE_SHARED(txtBuffer);
848 } else {
849 fputs(txtBuffer, fileHandle_p);
850 }
851 } /* if noError && writeCellData */
852
853 /************************************************************************/
854 /* write point data */
855
856 if (writePointData && Finley_noError()) {
857 /* mark the active data arrays */
858 bool_t set_scalar=FALSE, set_vector=FALSE, set_tensor=FALSE;
859 strcpy(txtBuffer, "<PointData");
860 for (dataIdx = 0; dataIdx<num_data; dataIdx++) {
861 if (!isEmpty(data_pp[dataIdx]) && !isCellCentered[dataIdx]) {
862 switch (getDataPointRank(data_pp[dataIdx])) {
863 case 0:
864 if (!set_scalar) {
865 strcat(txtBuffer, " Scalars=\"");
866 strcat(txtBuffer, names_p[dataIdx]);
867 strcat(txtBuffer, "\"");
868 set_scalar = TRUE;
869 }
870 break;
871 case 1:
872 if (!set_vector) {
873 strcat(txtBuffer, " Vectors=\"");
874 strcat(txtBuffer, names_p[dataIdx]);
875 strcat(txtBuffer, "\"");
876 set_vector = TRUE;
877 }
878 break;
879 case 2:
880 if (!set_tensor) {
881 strcat(txtBuffer, " Tensors=\"");
882 strcat(txtBuffer, names_p[dataIdx]);
883 strcat(txtBuffer, "\"");
884 set_tensor = TRUE;
885 }
886 break;
887 default:
888 sprintf(errorMsg, "saveVTK: data %s: VTK supports data with rank <= 2 only.", names_p[dataIdx]);
889 Finley_setError(VALUE_ERROR, errorMsg);
890 }
891 }
892 if (!Finley_noError())
893 break;
894 }
895 }
896 /* only continue if no error occurred */
897 if (writePointData && Finley_noError()) {
898 strcat(txtBuffer, ">\n");
899 if ( mpi_size > 1) {
900 MPI_RANK0_WRITE_SHARED(txtBuffer);
901 } else {
902 fputs(txtBuffer, fileHandle_p);
903 }
904
905 /* write the arrays */
906 for (dataIdx = 0; dataIdx < num_data; dataIdx++) {
907 if (!isEmpty(data_pp[dataIdx]) && !isCellCentered[dataIdx]) {
908 Finley_NodeMapping* nodeMapping;
909 dim_t rank = getDataPointRank(data_pp[dataIdx]);
910 dim_t nCompReqd = 1; /* number of components mpi_required */
911 if (getFunctionSpaceType(data_pp[dataIdx]) == FINLEY_REDUCED_NODES) {
912 nodeMapping = mesh_p->Nodes->reducedNodesMapping;
913 } else {
914 nodeMapping = mesh_p->Nodes->nodesMapping;
915 }
916 if (rank == 0) {
917 nCompReqd = 1;
918 shape = 0;
919 } else if (rank == 1) {
920 nCompReqd = 3;
921 shape = getDataPointShape(data_pp[dataIdx], 0);
922 if (shape > 3) {
923 Finley_setError(VALUE_ERROR, "saveVTK: rank 1 objects must have 3 components at most.");
924 }
925 } else {
926 nCompReqd = 9;
927 shape=getDataPointShape(data_pp[dataIdx], 0);
928 if (shape > 3 || shape != getDataPointShape(data_pp[dataIdx], 1)) {
929 Finley_setError(VALUE_ERROR, "saveVTK: rank 2 objects of shape 2x2 or 3x3 supported only.");
930 }
931 }
932 /* bail out if an error occurred */
933 if (!Finley_noError())
934 break;
935
936 sprintf(txtBuffer, tag_Float_DataArray, names_p[dataIdx], nCompReqd);
937 if ( mpi_size > 1) {
938 MPI_RANK0_WRITE_SHARED(txtBuffer);
939 } else {
940 fputs(txtBuffer, fileHandle_p);
941 }
942
943 txtBuffer[0] = '\0';
944 txtBufferInUse = 0;
945 for (i=0; i<mesh_p->Nodes->numNodes; i++) {
946 k = globalNodeIndex[i];
947 if ( (myFirstNode <= k) && (k < myLastNode) ) {
948 void* sampleBuffer=allocSampleBuffer(data_pp[dataIdx]);
949 __const double *values = getSampleDataRO(data_pp[dataIdx], nodeMapping->target[i], sampleBuffer);
950 /* if the number of mpi_required components is more than
951 * the number of actual components, pad with zeros.
952 * Probably only need to get shape of first element */
953 if (nCompReqd == 1) {
954 sprintf(tmpBuffer, SCALAR_FORMAT, values[0]);
955 } else if (nCompReqd == 3) {
956 if (shape==1) {
957 sprintf(tmpBuffer, VECTOR_FORMAT,
958 values[0], 0.f, 0.f);
959 } else if (shape==2) {
960 sprintf(tmpBuffer, VECTOR_FORMAT,
961 values[0], values[1], 0.f);
962 } else if (shape==3) {
963 sprintf(tmpBuffer, VECTOR_FORMAT,
964 values[0], values[1], values[2]);
965 }
966 } else if (nCompReqd == 9) {
967 if (shape==1) {
968 sprintf(tmpBuffer, TENSOR_FORMAT,
969 values[0], 0.f, 0.f,
970 0.f, 0.f, 0.f,
971 0.f, 0.f, 0.f);
972 } else if (shape==2) {
973 sprintf(tmpBuffer, TENSOR_FORMAT,
974 values[0], values[1], 0.f,
975 values[2], values[3], 0.f,
976 0.f, 0.f, 0.f);
977 } else if (shape==3) {
978 sprintf(tmpBuffer, TENSOR_FORMAT,
979 values[0], values[1], values[2],
980 values[3], values[4], values[5],
981 values[6], values[7], values[8]);
982 }
983 }
984 if ( mpi_size > 1) {
985 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
986 } else {
987 fputs(tmpBuffer, fileHandle_p);
988 }
989 freeSampleBuffer(sampleBuffer); /* no-one needs values anymore */
990 } /* if this is my node */
991 } /* for i (numNodes) */
992
993 if ( mpi_size > 1) {
994 MPI_WRITE_ORDERED(txtBuffer, txtBufferInUse);
995 MPI_RANK0_WRITE_SHARED(tag_End_DataArray);
996 } else {
997 fputs(tag_End_DataArray, fileHandle_p);
998 }
999 } /* !isEmpty && !isCellCentered */
1000 } /* for dataIdx */
1001
1002 strcpy(txtBuffer, "</PointData>\n");
1003 if ( mpi_size > 1) {
1004 MPI_RANK0_WRITE_SHARED(txtBuffer);
1005 } else {
1006 fputs(txtBuffer, fileHandle_p);
1007 }
1008 } /* if noError && writePointData */
1009
1010 /* Final write to VTK file */
1011 if (Finley_noError()) {
1012 if (mpi_size > 1) {
1013 MPI_RANK0_WRITE_SHARED(vtkFooter);
1014 } else {
1015 fputs(vtkFooter, fileHandle_p);
1016 }
1017 }
1018
1019 if ( mpi_size > 1) {
1020 #ifdef PASO_MPI
1021 MPI_File_close(&mpi_fileHandle_p);
1022 MPI_Barrier(mesh_p->Nodes->MPIInfo->comm);
1023 #endif
1024 } else {
1025 fclose(fileHandle_p);
1026 }
1027 TMPMEMFREE(isCellCentered);
1028 TMPMEMFREE(txtBuffer);
1029 }
1030

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26