/[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 2743 - (show annotations)
Mon Nov 16 05:28:47 2009 UTC (9 years, 9 months ago) by caltinay
File MIME type: text/plain
File size: 46122 byte(s)
Corrected memory allocation size for tensors in saveVTK. Fixes issue 460.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26