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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26