/[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 3145 - (show annotations)
Fri Sep 3 01:25:52 2010 UTC (8 years, 8 months ago) by jfenwick
File MIME type: text/plain
File size: 39349 byte(s)
Removed linearNodes and subElementNodes from ReferenceElementInfo
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14 /***************************************************************************/
15 /* Writes data and mesh in VTK XML format to a VTU file. */
16 /* Nodal data needs to be given on DUDLEY_NODES or DUDLEY_REDUCED_NODES */
17 /***************************************************************************/
18
19 #include "Mesh.h"
20 #include "Assemble.h"
21 #include "vtkCellType.h" /* copied from vtk source directory */
22 #include "paso/PasoUtil.h"
23
24 #define INT_FORMAT "%d "
25 #define LEN_INT_FORMAT (unsigned int)(9+2)
26 #define INT_NEWLINE_FORMAT "%d\n"
27 #define SCALAR_FORMAT "%12.6e\n"
28 #define VECTOR_FORMAT "%12.6e %12.6e %12.6e\n"
29 #define TENSOR_FORMAT "%12.6e %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e\n"
30 /* strlen("-1.234567e+789 ") == 15 */
31 #define LEN_TENSOR_FORMAT (unsigned int)(9*15+2)
32 #define NEWLINE "\n"
33 #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 return 1;
97 }
98
99 void Dudley_Mesh_saveVTK(const char *filename_p,
100 Dudley_Mesh *mesh_p,
101 const dim_t num_data,
102 char **names_p,
103 escriptDataC **data_pp,
104 const char* metadata,
105 const char*metadata_schema)
106 {
107 #ifdef PASO_MPI
108 MPI_File mpi_fileHandle_p;
109 MPI_Status mpi_status;
110 MPI_Request mpi_req;
111 MPI_Info mpi_info = MPI_INFO_NULL;
112 #endif
113 Paso_MPI_rank my_mpi_rank;
114 FILE *fileHandle_p = NULL;
115 char errorMsg[LenErrorMsg_MAX], *txtBuffer;
116 char tmpBuffer[LEN_TMP_BUFFER];
117 size_t txtBufferSize, txtBufferInUse, maxNameLen;
118 double *quadNodes_p = NULL;
119 dim_t dataIdx, nDim;
120 dim_t numCells=0, globalNumCells=0, numVTKNodesPerElement=0;
121 dim_t myNumPoints=0, globalNumPoints=0;
122 dim_t shape, NN=0, numCellFactor=1, myNumCells=0;
123 bool_t *isCellCentered;
124 bool_t writeCellData=FALSE, writePointData=FALSE, hasReducedElements=FALSE;
125 index_t myFirstNode=0, myLastNode=0, *globalNodeIndex=NULL;
126 index_t myFirstCell=0, k;
127 int mpi_size, i, j, l;
128 int cellType=0, nodeType=DUDLEY_NODES, elementType=DUDLEY_UNKNOWN;
129 Dudley_ElementFile *elements = NULL;
130 ElementTypeId typeId = NoRef;
131
132 const char *vtkHeader = \
133 "<?xml version=\"1.0\"?>\n" \
134 "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\"%s%s>\n%s%s" \
135 "<UnstructuredGrid>\n" \
136 "<Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n" \
137 "<Points>\n" \
138 "<DataArray NumberOfComponents=\"%d\" type=\"Float64\" format=\"ascii\">\n";
139 char *vtkFooter = "</Piece>\n</UnstructuredGrid>\n</VTKFile>\n";
140 const char *tag_Float_DataArray="<DataArray Name=\"%s\" type=\"Float64\" NumberOfComponents=\"%d\" format=\"ascii\">\n";
141 char *tags_End_Points_and_Start_Conn = "</DataArray>\n</Points>\n<Cells>\n<DataArray Name=\"connectivity\" type=\"Int32\" format=\"ascii\">\n" ;
142 char *tags_End_Conn_and_Start_Offset = "</DataArray>\n<DataArray Name=\"offsets\" type=\"Int32\" format=\"ascii\">\n";
143 char *tags_End_Offset_and_Start_Type = "</DataArray>\n<DataArray Name=\"types\" type=\"UInt8\" format=\"ascii\">\n";
144 char *tag_End_DataArray = "</DataArray>\n";
145
146 /* if there is no mesh we just return */
147 if (mesh_p==NULL) return;
148
149 nDim = mesh_p->Nodes->numDim;
150
151 if (nDim != 2 && nDim != 3) {
152 Dudley_setError(TYPE_ERROR, "saveVTK: spatial dimension 2 or 3 is supported only.");
153 return;
154 }
155 my_mpi_rank = mesh_p->Nodes->MPIInfo->rank;
156 mpi_size = mesh_p->Nodes->MPIInfo->size;
157
158 /************************************************************************
159 * open the file and check handle *
160 */
161 if (mpi_size > 1) {
162 #ifdef PASO_MPI
163 const int amode = MPI_MODE_CREATE|MPI_MODE_WRONLY|MPI_MODE_UNIQUE_OPEN;
164 int ierr;
165 if (my_mpi_rank == 0 && Paso_fileExists(filename_p)) {
166 remove(filename_p);
167 }
168 ierr = MPI_File_open(mesh_p->Nodes->MPIInfo->comm, (char*)filename_p,
169 amode, mpi_info, &mpi_fileHandle_p);
170 if (ierr != MPI_SUCCESS) {
171 sprintf(errorMsg, "saveVTK: File %s could not be opened for writing in parallel.", filename_p);
172 Dudley_setError(IO_ERROR, errorMsg);
173 } else {
174 ierr=MPI_File_set_view(mpi_fileHandle_p,MPI_DISPLACEMENT_CURRENT,
175 MPI_CHAR, MPI_CHAR, "native", mpi_info);
176 }
177 #endif /* PASO_MPI */
178 } else {
179 fileHandle_p = fopen(filename_p, "w");
180 if (fileHandle_p==NULL) {
181 sprintf(errorMsg, "saveVTK: File %s could not be opened for writing.", filename_p);
182 Dudley_setError(IO_ERROR, errorMsg);
183 }
184 }
185 if (!Paso_MPIInfo_noError(mesh_p->Nodes->MPIInfo)) return;
186
187 /* General note: From this point if an error occurs Dudley_setError is
188 * called and subsequent steps are skipped until the end of this function
189 * where allocated memory is freed and the file is closed. */
190
191 /************************************************************************/
192 /* find the mesh type to be written */
193
194 isCellCentered = TMPMEMALLOC(num_data, bool_t);
195 maxNameLen = 0;
196 if (!Dudley_checkPtr(isCellCentered)) {
197 for (dataIdx=0; dataIdx<num_data; ++dataIdx) {
198 if (! isEmpty(data_pp[dataIdx])) {
199 switch(getFunctionSpaceType(data_pp[dataIdx]) ) {
200 case DUDLEY_NODES:
201 nodeType = (nodeType == DUDLEY_REDUCED_NODES) ? DUDLEY_REDUCED_NODES : DUDLEY_NODES;
202 isCellCentered[dataIdx] = FALSE;
203 if (elementType==DUDLEY_UNKNOWN || elementType==DUDLEY_ELEMENTS) {
204 elementType = DUDLEY_ELEMENTS;
205 } else {
206 Dudley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
207 }
208 break;
209 case DUDLEY_REDUCED_NODES:
210 nodeType = DUDLEY_REDUCED_NODES;
211 isCellCentered[dataIdx] = FALSE;
212 if (elementType==DUDLEY_UNKNOWN || elementType==DUDLEY_ELEMENTS) {
213 elementType = DUDLEY_ELEMENTS;
214 } else {
215 Dudley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
216 }
217 break;
218 case DUDLEY_REDUCED_ELEMENTS:
219 hasReducedElements = TRUE;
220 case DUDLEY_ELEMENTS:
221 isCellCentered[dataIdx] = TRUE;
222 if (elementType==DUDLEY_UNKNOWN || elementType==DUDLEY_ELEMENTS) {
223 elementType = DUDLEY_ELEMENTS;
224 } else {
225 Dudley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
226 }
227 break;
228 case DUDLEY_REDUCED_FACE_ELEMENTS:
229 hasReducedElements = TRUE;
230 case DUDLEY_FACE_ELEMENTS:
231 isCellCentered[dataIdx] = TRUE;
232 if (elementType==DUDLEY_UNKNOWN || elementType==DUDLEY_FACE_ELEMENTS) {
233 elementType = DUDLEY_FACE_ELEMENTS;
234 } else {
235 Dudley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
236 }
237 break;
238 case DUDLEY_POINTS:
239 isCellCentered[dataIdx]=TRUE;
240 if (elementType==DUDLEY_UNKNOWN || elementType==DUDLEY_POINTS) {
241 elementType = DUDLEY_POINTS;
242 } else {
243 Dudley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
244 }
245 break;
246 default:
247 sprintf(errorMsg, "saveVTK: unknown function space type %d",getFunctionSpaceType(data_pp[dataIdx]));
248 Dudley_setError(TYPE_ERROR, errorMsg);
249 }
250 if (isCellCentered[dataIdx]) {
251 writeCellData = TRUE;
252 } else {
253 writePointData = TRUE;
254 }
255 maxNameLen = MAX(maxNameLen, strlen(names_p[dataIdx]));
256 }
257 }
258 }
259
260 /************************************************************************/
261 /* select number of points and the mesh component */
262
263 if (Dudley_noError()) {
264 if (nodeType == DUDLEY_REDUCED_NODES) {
265 myFirstNode = Dudley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
266 myLastNode = Dudley_NodeFile_getLastReducedNode(mesh_p->Nodes);
267 globalNumPoints = Dudley_NodeFile_getGlobalNumReducedNodes(mesh_p->Nodes);
268 globalNodeIndex = Dudley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
269 } else {
270 myFirstNode = Dudley_NodeFile_getFirstNode(mesh_p->Nodes);
271 myLastNode = Dudley_NodeFile_getLastNode(mesh_p->Nodes);
272 globalNumPoints = Dudley_NodeFile_getGlobalNumNodes(mesh_p->Nodes);
273 globalNodeIndex = Dudley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
274 }
275 myNumPoints = myLastNode - myFirstNode;
276 if (elementType==DUDLEY_UNKNOWN) elementType=DUDLEY_ELEMENTS;
277 switch(elementType) {
278 case DUDLEY_ELEMENTS:
279 elements = mesh_p->Elements;
280 break;
281 case DUDLEY_FACE_ELEMENTS:
282 elements = mesh_p->FaceElements;
283 break;
284 case DUDLEY_POINTS:
285 elements = mesh_p->Points;
286 break;
287 }
288 if (elements==NULL) {
289 Dudley_setError(SYSTEM_ERROR, "saveVTK: undefined element file");
290 } else {
291 /* map dudley element type to VTK element type */
292 numCells = elements->numElements;
293 globalNumCells = Dudley_ElementFile_getGlobalNumElements(elements);
294 myNumCells = Dudley_ElementFile_getMyNumElements(elements);
295 myFirstCell = Dudley_ElementFile_getFirstElement(elements);
296 NN = elements->numNodes;
297 if (hasReducedElements) {
298 quadNodes_p=elements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->QuadNodes;
299 } else {
300 quadNodes_p=elements->referenceElementSet->referenceElement->Parametrization->QuadNodes;
301 }
302 if (nodeType==DUDLEY_REDUCED_NODES) {
303 typeId = elements->referenceElementSet->referenceElement->LinearType->TypeId;
304 } else {
305 typeId = elements->referenceElementSet->referenceElement->Type->TypeId;
306 }
307 switch (typeId) {
308 case Point1:
309 case Line2Face:
310 cellType = VTK_VERTEX;
311 numVTKNodesPerElement = 1;
312 break;
313
314 case Line2:
315 case Tri3Face:
316 cellType = VTK_LINE;
317 numVTKNodesPerElement = 2;
318 break;
319
320 case Tri3:
321 case Tet4Face:
322 cellType = VTK_TRIANGLE;
323 numVTKNodesPerElement = 3;
324 break;
325
326 case Tet4:
327 cellType = VTK_TETRA;
328 numVTKNodesPerElement = 4;
329 break;
330
331 default:
332 sprintf(errorMsg, "saveVTK: Element type %s is not supported by VTK.", elements->referenceElementSet->referenceElement->Type->Name);
333 Dudley_setError(VALUE_ERROR, errorMsg);
334 }
335 }
336 }
337
338 /* allocate enough memory for text buffer */
339
340 txtBufferSize = strlen(vtkHeader) + 3*LEN_INT_FORMAT + (30+3*maxNameLen)+strlen(metadata)+strlen(metadata_schema);
341 if (mpi_size > 1) {
342 txtBufferSize = MAX(txtBufferSize, myNumPoints * LEN_TMP_BUFFER);
343 txtBufferSize = MAX(txtBufferSize, numCellFactor * myNumCells *
344 (LEN_INT_FORMAT * numVTKNodesPerElement + 1));
345 txtBufferSize = MAX(txtBufferSize,
346 numCellFactor * myNumCells * LEN_TENSOR_FORMAT);
347 txtBufferSize = MAX(txtBufferSize, myNumPoints * LEN_TENSOR_FORMAT);
348 }
349 txtBuffer = TMPMEMALLOC(txtBufferSize+1, char);
350
351 /* sets error if memory allocation failed */
352 Dudley_checkPtr(txtBuffer);
353
354 /************************************************************************/
355 /* write number of points and the mesh component */
356
357 if (Dudley_noError()) {
358 int flag1=0;
359 if (DUDLEY_REDUCED_NODES == nodeType) {
360 flag1=1;
361 } else if (numVTKNodesPerElement != elements->referenceElementSet->referenceElement->Type->numNodes) {
362 flag1=1;
363 }
364 if (strlen(metadata)>0) {
365 if (strlen(metadata_schema)>0) {
366 sprintf(txtBuffer, vtkHeader," ",metadata_schema,metadata,"\n",globalNumPoints, numCellFactor*globalNumCells, 3);
367 } else {
368 sprintf(txtBuffer, vtkHeader,"","",metadata,"\n",globalNumPoints, numCellFactor*globalNumCells, 3);
369 }
370 } else {
371 if (strlen(metadata_schema)>0) {
372 sprintf(txtBuffer, vtkHeader," ",metadata_schema,"","",globalNumPoints, numCellFactor*globalNumCells, 3);
373 } else {
374 sprintf(txtBuffer, vtkHeader,"","","","",globalNumPoints, numCellFactor*globalNumCells, 3);
375 }
376 }
377
378 if (mpi_size > 1) {
379 /* write the nodes */
380 MPI_RANK0_WRITE_SHARED(txtBuffer);
381 txtBuffer[0] = '\0';
382 txtBufferInUse = 0;
383 if (nDim==2) {
384 for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
385 if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
386 sprintf(tmpBuffer, VECTOR_FORMAT,
387 mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
388 mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
389 0.);
390 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
391 }
392 }
393 } else {
394 for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
395 if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
396 sprintf(tmpBuffer, VECTOR_FORMAT,
397 mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
398 mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
399 mesh_p->Nodes->Coordinates[INDEX2(2, i, nDim)]);
400 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
401 }
402 }
403 } /* nDim */
404 MPI_WRITE_ORDERED(txtBuffer);
405
406 /* write the cells */
407 MPI_RANK0_WRITE_SHARED(tags_End_Points_and_Start_Conn);
408 txtBuffer[0] = '\0';
409 txtBufferInUse = 0;
410 if (!flag1) {
411 for (i = 0; i < numCells; i++) {
412 if (elements->Owner[i] == my_mpi_rank) {
413 for (j = 0; j < numVTKNodesPerElement; j++) {
414 sprintf(tmpBuffer, INT_FORMAT, globalNodeIndex[elements->Nodes[INDEX2(j, i, NN)]]);
415 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
416 }
417 __STRCAT(txtBuffer, NEWLINE, txtBufferInUse);
418 }
419 }
420 } else {
421 for (i = 0; i < numCells; i++) {
422 if (elements->Owner[i] == my_mpi_rank) {
423 for (l = 0; l < numCellFactor; l++) {
424 const int idx=l*numVTKNodesPerElement;
425 for (j = 0; j < numVTKNodesPerElement; j++) {
426 sprintf(tmpBuffer, INT_FORMAT, globalNodeIndex[elements->Nodes[INDEX2(idx+j, i, NN)]]);
427 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
428 }
429 __STRCAT(txtBuffer, NEWLINE, txtBufferInUse);
430 }
431 }
432 }
433 } /* nodeIndex */
434 MPI_WRITE_ORDERED(txtBuffer);
435
436 /* write the offsets */
437 MPI_RANK0_WRITE_SHARED(tags_End_Conn_and_Start_Offset);
438 txtBuffer[0] = '\0';
439 txtBufferInUse = 0;
440 for (i = numVTKNodesPerElement*(myFirstCell*numCellFactor+1);
441 i <= (myFirstCell+myNumCells)*numVTKNodesPerElement*numCellFactor; i += numVTKNodesPerElement)
442 {
443 sprintf(tmpBuffer, INT_NEWLINE_FORMAT, i);
444 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
445 }
446 MPI_WRITE_ORDERED(txtBuffer);
447
448 /* write element type */
449 sprintf(tmpBuffer, INT_NEWLINE_FORMAT, cellType);
450 MPI_RANK0_WRITE_SHARED(tags_End_Offset_and_Start_Type);
451 txtBuffer[0] = '\0';
452 txtBufferInUse = 0;
453 for (i = numVTKNodesPerElement*(myFirstCell*numCellFactor+1);
454 i <= (myFirstCell+myNumCells)*numVTKNodesPerElement*numCellFactor; i += numVTKNodesPerElement)
455 {
456 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
457 }
458 MPI_WRITE_ORDERED(txtBuffer);
459 /* finalize cell information */
460 strcpy(txtBuffer, "</DataArray>\n</Cells>\n");
461 MPI_RANK0_WRITE_SHARED(txtBuffer);
462
463 } else { /***** mpi_size == 1 *****/
464
465 /* write the nodes */
466 fputs(txtBuffer, fileHandle_p);
467 if (nDim==2) {
468 for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
469 if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
470 fprintf(fileHandle_p, VECTOR_FORMAT,
471 mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
472 mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
473 0.);
474 }
475 }
476 } else {
477 for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
478 if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
479 fprintf(fileHandle_p, VECTOR_FORMAT,
480 mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
481 mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
482 mesh_p->Nodes->Coordinates[INDEX2(2, i, nDim)]);
483 }
484 }
485 } /* nDim */
486
487 /* write the cells */
488 fputs(tags_End_Points_and_Start_Conn, fileHandle_p);
489 if (!flag1) {
490 for (i = 0; i < numCells; i++) {
491 for (j = 0; j < numVTKNodesPerElement; j++) {
492 fprintf(fileHandle_p, INT_FORMAT, globalNodeIndex[elements->Nodes[INDEX2(j, i, NN)]]);
493 }
494 fprintf(fileHandle_p, NEWLINE);
495 }
496 } else {
497 for (i = 0; i < numCells; i++) {
498 for (l = 0; l < numCellFactor; l++) {
499 const int idx=l*numVTKNodesPerElement;
500 for (j = 0; j < numVTKNodesPerElement; j++) {
501 fprintf(fileHandle_p, INT_FORMAT, globalNodeIndex[elements->Nodes[INDEX2(idx+j, i, NN)]]);
502 }
503 fprintf(fileHandle_p, NEWLINE);
504 }
505 }
506 } /* nodeIndex */
507
508 /* write the offsets */
509 fputs(tags_End_Conn_and_Start_Offset, fileHandle_p);
510 for (i = numVTKNodesPerElement; i <= numCells*numVTKNodesPerElement*numCellFactor; i += numVTKNodesPerElement) {
511 fprintf(fileHandle_p, INT_NEWLINE_FORMAT, i);
512 }
513
514 /* write element type */
515 sprintf(tmpBuffer, INT_NEWLINE_FORMAT, cellType);
516 fputs(tags_End_Offset_and_Start_Type, fileHandle_p);
517 for (i = 0; i < numCells*numCellFactor; i++)
518 fputs(tmpBuffer, fileHandle_p);
519 /* finalize cell information */
520 fputs("</DataArray>\n</Cells>\n", fileHandle_p);
521 } /* mpi_size */
522
523 } /* Dudley_noError */
524
525 /************************************************************************/
526 /* write cell data */
527
528 if (writeCellData && Dudley_noError()) {
529 bool_t set_scalar=FALSE, set_vector=FALSE, set_tensor=FALSE;
530 /* mark the active data arrays */
531 strcpy(txtBuffer, "<CellData");
532 for (dataIdx = 0; dataIdx < num_data; dataIdx++) {
533 if (!isEmpty(data_pp[dataIdx]) && isCellCentered[dataIdx]) {
534 /* rank == 0 <--> scalar data */
535 /* rank == 1 <--> vector data */
536 /* rank == 2 <--> tensor data */
537 switch (getDataPointRank(data_pp[dataIdx])) {
538 case 0:
539 if (!set_scalar) {
540 strcat(txtBuffer, " Scalars=\"");
541 strcat(txtBuffer, names_p[dataIdx]);
542 strcat(txtBuffer, "\"");
543 set_scalar = TRUE;
544 }
545 break;
546 case 1:
547 if (!set_vector) {
548 strcat(txtBuffer, " Vectors=\"");
549 strcat(txtBuffer, names_p[dataIdx]);
550 strcat(txtBuffer, "\"");
551 set_vector = TRUE;
552 }
553 break;
554 case 2:
555 if (!set_tensor) {
556 strcat(txtBuffer, " Tensors=\"");
557 strcat(txtBuffer, names_p[dataIdx]);
558 strcat(txtBuffer, "\"");
559 set_tensor = TRUE;
560 }
561 break;
562 default:
563 sprintf(errorMsg, "saveVTK: data %s: VTK supports data with rank <= 2 only.", names_p[dataIdx]);
564 Dudley_setError(VALUE_ERROR, errorMsg);
565 }
566 }
567 if (!Dudley_noError())
568 break;
569 }
570 }
571 /* only continue if no error occurred */
572 if (writeCellData && Dudley_noError()) {
573 strcat(txtBuffer, ">\n");
574 if ( mpi_size > 1) {
575 MPI_RANK0_WRITE_SHARED(txtBuffer);
576 } else {
577 fputs(txtBuffer, fileHandle_p);
578 }
579
580 /* write the arrays */
581 for (dataIdx = 0; dataIdx < num_data; dataIdx++) {
582 if (!isEmpty(data_pp[dataIdx]) && isCellCentered[dataIdx]) {
583 dim_t numPointsPerSample=getNumDataPointsPerSample(data_pp[dataIdx]);
584 dim_t rank = getDataPointRank(data_pp[dataIdx]);
585 dim_t nComp = getDataPointSize(data_pp[dataIdx]);
586 dim_t nCompReqd = 1; /* number of components mpi_required */
587 if (rank == 0) {
588 nCompReqd = 1;
589 shape = 0;
590 } else if (rank == 1) {
591 nCompReqd = 3;
592 shape = getDataPointShape(data_pp[dataIdx], 0);
593 if (shape > 3) {
594 Dudley_setError(VALUE_ERROR, "saveVTK: rank 1 objects must have 3 components at most.");
595 }
596 } else {
597 nCompReqd = 9;
598 shape = getDataPointShape(data_pp[dataIdx], 0);
599 if (shape > 3 || shape != getDataPointShape(data_pp[dataIdx], 1)) {
600 Dudley_setError(VALUE_ERROR, "saveVTK: rank 2 objects of shape 2x2 or 3x3 supported only.");
601 }
602 }
603 /* bail out if an error occurred */
604 if (!Dudley_noError())
605 break;
606
607 sprintf(txtBuffer, tag_Float_DataArray, names_p[dataIdx], nCompReqd);
608 if ( mpi_size > 1) {
609 MPI_RANK0_WRITE_SHARED(txtBuffer);
610 } else {
611 fputs(txtBuffer, fileHandle_p);
612 }
613
614 txtBuffer[0] = '\0';
615 txtBufferInUse = 0;
616 for (i=0; i<numCells; i++) {
617 if (elements->Owner[i] == my_mpi_rank) {
618 __const double *values = getSampleDataRO(data_pp[dataIdx], i);
619 for (l = 0; l < numCellFactor; l++) {
620 double sampleAvg[NCOMP_MAX];
621 dim_t nCompUsed = MIN(nComp, NCOMP_MAX);
622
623 /* average over number of points in the sample */
624 if (isExpanded(data_pp[dataIdx])) {
625 dim_t hits=0;
626 for (k=0; k<nCompUsed; k++) sampleAvg[k]=0;
627 for (j=0; j<numPointsPerSample; j++) {
628 if (nodeInQuadrant(quadNodes_p, typeId, j, l)) {
629 hits++;
630 for (k=0; k<nCompUsed; k++) {
631 sampleAvg[k] += values[INDEX2(k,j,nComp)];
632 }
633 }
634 }
635 for (k=0; k<nCompUsed; k++)
636 sampleAvg[k] /= MAX(hits, 1);
637 } else {
638 for (k=0; k<nCompUsed; k++)
639 sampleAvg[k] = values[k];
640 } /* isExpanded */
641
642 /* if the number of required components is more than
643 * the number of actual components, pad with zeros
644 */
645 /* probably only need to get shape of first element */
646 if (nCompReqd == 1) {
647 sprintf(tmpBuffer, SCALAR_FORMAT, sampleAvg[0]);
648 } else if (nCompReqd == 3) {
649 if (shape==1) {
650 sprintf(tmpBuffer, VECTOR_FORMAT,
651 sampleAvg[0], 0.f, 0.f);
652 } else if (shape==2) {
653 sprintf(tmpBuffer, VECTOR_FORMAT,
654 sampleAvg[0], sampleAvg[1], 0.f);
655 } else if (shape==3) {
656 sprintf(tmpBuffer, VECTOR_FORMAT,
657 sampleAvg[0],sampleAvg[1],sampleAvg[2]);
658 }
659 } else if (nCompReqd == 9) {
660 if (shape==1) {
661 sprintf(tmpBuffer, TENSOR_FORMAT,
662 sampleAvg[0], 0.f, 0.f,
663 0.f, 0.f, 0.f,
664 0.f, 0.f, 0.f);
665 } else if (shape==2) {
666 sprintf(tmpBuffer, TENSOR_FORMAT,
667 sampleAvg[0], sampleAvg[1], 0.f,
668 sampleAvg[2], sampleAvg[3], 0.f,
669 0.f, 0.f, 0.f);
670 } else if (shape==3) {
671 sprintf(tmpBuffer, TENSOR_FORMAT,
672 sampleAvg[0],sampleAvg[1],sampleAvg[2],
673 sampleAvg[3],sampleAvg[4],sampleAvg[5],
674 sampleAvg[6],sampleAvg[7],sampleAvg[8]);
675 }
676 }
677 if ( mpi_size > 1) {
678 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
679 } else {
680 fputs(tmpBuffer, fileHandle_p);
681 }
682 } /* for l (numCellFactor) */
683 } /* if I am the owner */
684 } /* for i (numCells) */
685
686 if ( mpi_size > 1) {
687 MPI_WRITE_ORDERED(txtBuffer);
688 MPI_RANK0_WRITE_SHARED(tag_End_DataArray);
689 } else {
690 fputs(tag_End_DataArray, fileHandle_p);
691 }
692 } /* !isEmpty && cellCentered */
693 } /* for dataIdx */
694
695 strcpy(txtBuffer, "</CellData>\n");
696 if ( mpi_size > 1) {
697 MPI_RANK0_WRITE_SHARED(txtBuffer);
698 } else {
699 fputs(txtBuffer, fileHandle_p);
700 }
701 } /* if noError && writeCellData */
702
703 /************************************************************************/
704 /* write point data */
705
706 if (writePointData && Dudley_noError()) {
707 /* mark the active data arrays */
708 bool_t set_scalar=FALSE, set_vector=FALSE, set_tensor=FALSE;
709 strcpy(txtBuffer, "<PointData");
710 for (dataIdx = 0; dataIdx<num_data; dataIdx++) {
711 if (!isEmpty(data_pp[dataIdx]) && !isCellCentered[dataIdx]) {
712 switch (getDataPointRank(data_pp[dataIdx])) {
713 case 0:
714 if (!set_scalar) {
715 strcat(txtBuffer, " Scalars=\"");
716 strcat(txtBuffer, names_p[dataIdx]);
717 strcat(txtBuffer, "\"");
718 set_scalar = TRUE;
719 }
720 break;
721 case 1:
722 if (!set_vector) {
723 strcat(txtBuffer, " Vectors=\"");
724 strcat(txtBuffer, names_p[dataIdx]);
725 strcat(txtBuffer, "\"");
726 set_vector = TRUE;
727 }
728 break;
729 case 2:
730 if (!set_tensor) {
731 strcat(txtBuffer, " Tensors=\"");
732 strcat(txtBuffer, names_p[dataIdx]);
733 strcat(txtBuffer, "\"");
734 set_tensor = TRUE;
735 }
736 break;
737 default:
738 sprintf(errorMsg, "saveVTK: data %s: VTK supports data with rank <= 2 only.", names_p[dataIdx]);
739 Dudley_setError(VALUE_ERROR, errorMsg);
740 }
741 }
742 if (!Dudley_noError())
743 break;
744 }
745 }
746 /* only continue if no error occurred */
747 if (writePointData && Dudley_noError()) {
748 strcat(txtBuffer, ">\n");
749 if ( mpi_size > 1) {
750 MPI_RANK0_WRITE_SHARED(txtBuffer);
751 } else {
752 fputs(txtBuffer, fileHandle_p);
753 }
754
755 /* write the arrays */
756 for (dataIdx = 0; dataIdx < num_data; dataIdx++) {
757 if (!isEmpty(data_pp[dataIdx]) && !isCellCentered[dataIdx]) {
758 Dudley_NodeMapping* nodeMapping;
759 dim_t rank = getDataPointRank(data_pp[dataIdx]);
760 dim_t nCompReqd = 1; /* number of components mpi_required */
761 if (getFunctionSpaceType(data_pp[dataIdx]) == DUDLEY_REDUCED_NODES) {
762 nodeMapping = mesh_p->Nodes->reducedNodesMapping;
763 } else {
764 nodeMapping = mesh_p->Nodes->nodesMapping;
765 }
766 if (rank == 0) {
767 nCompReqd = 1;
768 shape = 0;
769 } else if (rank == 1) {
770 nCompReqd = 3;
771 shape = getDataPointShape(data_pp[dataIdx], 0);
772 if (shape > 3) {
773 Dudley_setError(VALUE_ERROR, "saveVTK: rank 1 objects must have 3 components at most.");
774 }
775 } else {
776 nCompReqd = 9;
777 shape=getDataPointShape(data_pp[dataIdx], 0);
778 if (shape > 3 || shape != getDataPointShape(data_pp[dataIdx], 1)) {
779 Dudley_setError(VALUE_ERROR, "saveVTK: rank 2 objects of shape 2x2 or 3x3 supported only.");
780 }
781 }
782 /* bail out if an error occurred */
783 if (!Dudley_noError())
784 break;
785
786 sprintf(txtBuffer, tag_Float_DataArray, names_p[dataIdx], nCompReqd);
787 if ( mpi_size > 1) {
788 MPI_RANK0_WRITE_SHARED(txtBuffer);
789 } else {
790 fputs(txtBuffer, fileHandle_p);
791 }
792
793 txtBuffer[0] = '\0';
794 txtBufferInUse = 0;
795 for (i=0; i<mesh_p->Nodes->numNodes; i++) {
796 k = globalNodeIndex[i];
797 if ( (myFirstNode <= k) && (k < myLastNode) ) {
798 __const double *values = getSampleDataRO(data_pp[dataIdx], nodeMapping->target[i]);
799 /* if the number of mpi_required components is more than
800 * the number of actual components, pad with zeros.
801 * Probably only need to get shape of first element */
802 if (nCompReqd == 1) {
803 sprintf(tmpBuffer, SCALAR_FORMAT, values[0]);
804 } else if (nCompReqd == 3) {
805 if (shape==1) {
806 sprintf(tmpBuffer, VECTOR_FORMAT,
807 values[0], 0.f, 0.f);
808 } else if (shape==2) {
809 sprintf(tmpBuffer, VECTOR_FORMAT,
810 values[0], values[1], 0.f);
811 } else if (shape==3) {
812 sprintf(tmpBuffer, VECTOR_FORMAT,
813 values[0], values[1], values[2]);
814 }
815 } else if (nCompReqd == 9) {
816 if (shape==1) {
817 sprintf(tmpBuffer, TENSOR_FORMAT,
818 values[0], 0.f, 0.f,
819 0.f, 0.f, 0.f,
820 0.f, 0.f, 0.f);
821 } else if (shape==2) {
822 sprintf(tmpBuffer, TENSOR_FORMAT,
823 values[0], values[1], 0.f,
824 values[2], values[3], 0.f,
825 0.f, 0.f, 0.f);
826 } else if (shape==3) {
827 sprintf(tmpBuffer, TENSOR_FORMAT,
828 values[0], values[1], values[2],
829 values[3], values[4], values[5],
830 values[6], values[7], values[8]);
831 }
832 }
833 if ( mpi_size > 1) {
834 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
835 } else {
836 fputs(tmpBuffer, fileHandle_p);
837 }
838 } /* if this is my node */
839 } /* for i (numNodes) */
840
841 if ( mpi_size > 1) {
842 MPI_WRITE_ORDERED(txtBuffer);
843 MPI_RANK0_WRITE_SHARED(tag_End_DataArray);
844 } else {
845 fputs(tag_End_DataArray, fileHandle_p);
846 }
847 } /* !isEmpty && !isCellCentered */
848 } /* for dataIdx */
849
850 strcpy(txtBuffer, "</PointData>\n");
851 if ( mpi_size > 1) {
852 MPI_RANK0_WRITE_SHARED(txtBuffer);
853 } else {
854 fputs(txtBuffer, fileHandle_p);
855 }
856 } /* if noError && writePointData */
857
858 /* Final write to VTK file */
859 if (Dudley_noError()) {
860 if (mpi_size > 1) {
861 MPI_RANK0_WRITE_SHARED(vtkFooter);
862 } else {
863 fputs(vtkFooter, fileHandle_p);
864 }
865 }
866
867 if ( mpi_size > 1) {
868 #ifdef PASO_MPI
869 MPI_File_close(&mpi_fileHandle_p);
870 MPI_Barrier(mesh_p->Nodes->MPIInfo->comm);
871 #endif
872 } else {
873 fclose(fileHandle_p);
874 }
875 TMPMEMFREE(isCellCentered);
876 TMPMEMFREE(txtBuffer);
877 }
878

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26