/[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 3126 - (show annotations)
Wed Sep 1 00:37:53 2010 UTC (8 years, 8 months ago) by jfenwick
File MIME type: text/plain
File size: 39634 byte(s)
Updated referenceElements and enums
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 const index_t *nodeIndex;
359 if (DUDLEY_REDUCED_NODES == nodeType) {
360 nodeIndex = elements->referenceElementSet->referenceElement->Type->linearNodes;
361 } else if (numVTKNodesPerElement != elements->referenceElementSet->referenceElement->Type->numNodes) {
362 nodeIndex = elements->referenceElementSet->referenceElement->Type->relevantGeoNodes;
363 } else {
364 nodeIndex = NULL;
365 }
366
367 if (strlen(metadata)>0) {
368 if (strlen(metadata_schema)>0) {
369 sprintf(txtBuffer, vtkHeader," ",metadata_schema,metadata,"\n",globalNumPoints, numCellFactor*globalNumCells, 3);
370 } else {
371 sprintf(txtBuffer, vtkHeader,"","",metadata,"\n",globalNumPoints, numCellFactor*globalNumCells, 3);
372 }
373 } else {
374 if (strlen(metadata_schema)>0) {
375 sprintf(txtBuffer, vtkHeader," ",metadata_schema,"","",globalNumPoints, numCellFactor*globalNumCells, 3);
376 } else {
377 sprintf(txtBuffer, vtkHeader,"","","","",globalNumPoints, numCellFactor*globalNumCells, 3);
378 }
379 }
380
381 if (mpi_size > 1) {
382 /* write the nodes */
383 MPI_RANK0_WRITE_SHARED(txtBuffer);
384 txtBuffer[0] = '\0';
385 txtBufferInUse = 0;
386 if (nDim==2) {
387 for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
388 if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
389 sprintf(tmpBuffer, VECTOR_FORMAT,
390 mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
391 mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
392 0.);
393 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
394 }
395 }
396 } else {
397 for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
398 if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
399 sprintf(tmpBuffer, VECTOR_FORMAT,
400 mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
401 mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
402 mesh_p->Nodes->Coordinates[INDEX2(2, i, nDim)]);
403 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
404 }
405 }
406 } /* nDim */
407 MPI_WRITE_ORDERED(txtBuffer);
408
409 /* write the cells */
410 MPI_RANK0_WRITE_SHARED(tags_End_Points_and_Start_Conn);
411 txtBuffer[0] = '\0';
412 txtBufferInUse = 0;
413 if (nodeIndex == NULL) {
414 for (i = 0; i < numCells; i++) {
415 if (elements->Owner[i] == my_mpi_rank) {
416 for (j = 0; j < numVTKNodesPerElement; j++) {
417 sprintf(tmpBuffer, INT_FORMAT, globalNodeIndex[elements->Nodes[INDEX2(j, i, NN)]]);
418 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
419 }
420 __STRCAT(txtBuffer, NEWLINE, txtBufferInUse);
421 }
422 }
423 } else {
424 for (i = 0; i < numCells; i++) {
425 if (elements->Owner[i] == my_mpi_rank) {
426 for (l = 0; l < numCellFactor; l++) {
427 const int* idx=&nodeIndex[l*numVTKNodesPerElement];
428 for (j = 0; j < numVTKNodesPerElement; j++) {
429 sprintf(tmpBuffer, INT_FORMAT, globalNodeIndex[elements->Nodes[INDEX2(idx[j], i, NN)]]);
430 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
431 }
432 __STRCAT(txtBuffer, NEWLINE, txtBufferInUse);
433 }
434 }
435 }
436 } /* nodeIndex */
437 MPI_WRITE_ORDERED(txtBuffer);
438
439 /* write the offsets */
440 MPI_RANK0_WRITE_SHARED(tags_End_Conn_and_Start_Offset);
441 txtBuffer[0] = '\0';
442 txtBufferInUse = 0;
443 for (i = numVTKNodesPerElement*(myFirstCell*numCellFactor+1);
444 i <= (myFirstCell+myNumCells)*numVTKNodesPerElement*numCellFactor; i += numVTKNodesPerElement)
445 {
446 sprintf(tmpBuffer, INT_NEWLINE_FORMAT, i);
447 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
448 }
449 MPI_WRITE_ORDERED(txtBuffer);
450
451 /* write element type */
452 sprintf(tmpBuffer, INT_NEWLINE_FORMAT, cellType);
453 MPI_RANK0_WRITE_SHARED(tags_End_Offset_and_Start_Type);
454 txtBuffer[0] = '\0';
455 txtBufferInUse = 0;
456 for (i = numVTKNodesPerElement*(myFirstCell*numCellFactor+1);
457 i <= (myFirstCell+myNumCells)*numVTKNodesPerElement*numCellFactor; i += numVTKNodesPerElement)
458 {
459 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
460 }
461 MPI_WRITE_ORDERED(txtBuffer);
462 /* finalize cell information */
463 strcpy(txtBuffer, "</DataArray>\n</Cells>\n");
464 MPI_RANK0_WRITE_SHARED(txtBuffer);
465
466 } else { /***** mpi_size == 1 *****/
467
468 /* write the nodes */
469 fputs(txtBuffer, fileHandle_p);
470 if (nDim==2) {
471 for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
472 if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
473 fprintf(fileHandle_p, VECTOR_FORMAT,
474 mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
475 mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
476 0.);
477 }
478 }
479 } else {
480 for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
481 if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
482 fprintf(fileHandle_p, VECTOR_FORMAT,
483 mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
484 mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
485 mesh_p->Nodes->Coordinates[INDEX2(2, i, nDim)]);
486 }
487 }
488 } /* nDim */
489
490 /* write the cells */
491 fputs(tags_End_Points_and_Start_Conn, fileHandle_p);
492 if (nodeIndex == NULL) {
493 for (i = 0; i < numCells; i++) {
494 for (j = 0; j < numVTKNodesPerElement; j++) {
495 fprintf(fileHandle_p, INT_FORMAT, globalNodeIndex[elements->Nodes[INDEX2(j, i, NN)]]);
496 }
497 fprintf(fileHandle_p, NEWLINE);
498 }
499 } else {
500 for (i = 0; i < numCells; i++) {
501 for (l = 0; l < numCellFactor; l++) {
502 const int* idx=&nodeIndex[l*numVTKNodesPerElement];
503 for (j = 0; j < numVTKNodesPerElement; j++) {
504 fprintf(fileHandle_p, INT_FORMAT, globalNodeIndex[elements->Nodes[INDEX2(idx[j], i, NN)]]);
505 }
506 fprintf(fileHandle_p, NEWLINE);
507 }
508 }
509 } /* nodeIndex */
510
511 /* write the offsets */
512 fputs(tags_End_Conn_and_Start_Offset, fileHandle_p);
513 for (i = numVTKNodesPerElement; i <= numCells*numVTKNodesPerElement*numCellFactor; i += numVTKNodesPerElement) {
514 fprintf(fileHandle_p, INT_NEWLINE_FORMAT, i);
515 }
516
517 /* write element type */
518 sprintf(tmpBuffer, INT_NEWLINE_FORMAT, cellType);
519 fputs(tags_End_Offset_and_Start_Type, fileHandle_p);
520 for (i = 0; i < numCells*numCellFactor; i++)
521 fputs(tmpBuffer, fileHandle_p);
522 /* finalize cell information */
523 fputs("</DataArray>\n</Cells>\n", fileHandle_p);
524 } /* mpi_size */
525
526 } /* Dudley_noError */
527
528 /************************************************************************/
529 /* write cell data */
530
531 if (writeCellData && Dudley_noError()) {
532 bool_t set_scalar=FALSE, set_vector=FALSE, set_tensor=FALSE;
533 /* mark the active data arrays */
534 strcpy(txtBuffer, "<CellData");
535 for (dataIdx = 0; dataIdx < num_data; dataIdx++) {
536 if (!isEmpty(data_pp[dataIdx]) && isCellCentered[dataIdx]) {
537 /* rank == 0 <--> scalar data */
538 /* rank == 1 <--> vector data */
539 /* rank == 2 <--> tensor data */
540 switch (getDataPointRank(data_pp[dataIdx])) {
541 case 0:
542 if (!set_scalar) {
543 strcat(txtBuffer, " Scalars=\"");
544 strcat(txtBuffer, names_p[dataIdx]);
545 strcat(txtBuffer, "\"");
546 set_scalar = TRUE;
547 }
548 break;
549 case 1:
550 if (!set_vector) {
551 strcat(txtBuffer, " Vectors=\"");
552 strcat(txtBuffer, names_p[dataIdx]);
553 strcat(txtBuffer, "\"");
554 set_vector = TRUE;
555 }
556 break;
557 case 2:
558 if (!set_tensor) {
559 strcat(txtBuffer, " Tensors=\"");
560 strcat(txtBuffer, names_p[dataIdx]);
561 strcat(txtBuffer, "\"");
562 set_tensor = TRUE;
563 }
564 break;
565 default:
566 sprintf(errorMsg, "saveVTK: data %s: VTK supports data with rank <= 2 only.", names_p[dataIdx]);
567 Dudley_setError(VALUE_ERROR, errorMsg);
568 }
569 }
570 if (!Dudley_noError())
571 break;
572 }
573 }
574 /* only continue if no error occurred */
575 if (writeCellData && Dudley_noError()) {
576 strcat(txtBuffer, ">\n");
577 if ( mpi_size > 1) {
578 MPI_RANK0_WRITE_SHARED(txtBuffer);
579 } else {
580 fputs(txtBuffer, fileHandle_p);
581 }
582
583 /* write the arrays */
584 for (dataIdx = 0; dataIdx < num_data; dataIdx++) {
585 if (!isEmpty(data_pp[dataIdx]) && isCellCentered[dataIdx]) {
586 dim_t numPointsPerSample=getNumDataPointsPerSample(data_pp[dataIdx]);
587 dim_t rank = getDataPointRank(data_pp[dataIdx]);
588 dim_t nComp = getDataPointSize(data_pp[dataIdx]);
589 dim_t nCompReqd = 1; /* number of components mpi_required */
590 if (rank == 0) {
591 nCompReqd = 1;
592 shape = 0;
593 } else if (rank == 1) {
594 nCompReqd = 3;
595 shape = getDataPointShape(data_pp[dataIdx], 0);
596 if (shape > 3) {
597 Dudley_setError(VALUE_ERROR, "saveVTK: rank 1 objects must have 3 components at most.");
598 }
599 } else {
600 nCompReqd = 9;
601 shape = getDataPointShape(data_pp[dataIdx], 0);
602 if (shape > 3 || shape != getDataPointShape(data_pp[dataIdx], 1)) {
603 Dudley_setError(VALUE_ERROR, "saveVTK: rank 2 objects of shape 2x2 or 3x3 supported only.");
604 }
605 }
606 /* bail out if an error occurred */
607 if (!Dudley_noError())
608 break;
609
610 sprintf(txtBuffer, tag_Float_DataArray, names_p[dataIdx], nCompReqd);
611 if ( mpi_size > 1) {
612 MPI_RANK0_WRITE_SHARED(txtBuffer);
613 } else {
614 fputs(txtBuffer, fileHandle_p);
615 }
616
617 txtBuffer[0] = '\0';
618 txtBufferInUse = 0;
619 for (i=0; i<numCells; i++) {
620 if (elements->Owner[i] == my_mpi_rank) {
621 __const double *values = getSampleDataRO(data_pp[dataIdx], i);
622 for (l = 0; l < numCellFactor; l++) {
623 double sampleAvg[NCOMP_MAX];
624 dim_t nCompUsed = MIN(nComp, NCOMP_MAX);
625
626 /* average over number of points in the sample */
627 if (isExpanded(data_pp[dataIdx])) {
628 dim_t hits=0;
629 for (k=0; k<nCompUsed; k++) sampleAvg[k]=0;
630 for (j=0; j<numPointsPerSample; j++) {
631 if (nodeInQuadrant(quadNodes_p, typeId, j, l)) {
632 hits++;
633 for (k=0; k<nCompUsed; k++) {
634 sampleAvg[k] += values[INDEX2(k,j,nComp)];
635 }
636 }
637 }
638 for (k=0; k<nCompUsed; k++)
639 sampleAvg[k] /= MAX(hits, 1);
640 } else {
641 for (k=0; k<nCompUsed; k++)
642 sampleAvg[k] = values[k];
643 } /* isExpanded */
644
645 /* if the number of required components is more than
646 * the number of actual components, pad with zeros
647 */
648 /* probably only need to get shape of first element */
649 if (nCompReqd == 1) {
650 sprintf(tmpBuffer, SCALAR_FORMAT, sampleAvg[0]);
651 } else if (nCompReqd == 3) {
652 if (shape==1) {
653 sprintf(tmpBuffer, VECTOR_FORMAT,
654 sampleAvg[0], 0.f, 0.f);
655 } else if (shape==2) {
656 sprintf(tmpBuffer, VECTOR_FORMAT,
657 sampleAvg[0], sampleAvg[1], 0.f);
658 } else if (shape==3) {
659 sprintf(tmpBuffer, VECTOR_FORMAT,
660 sampleAvg[0],sampleAvg[1],sampleAvg[2]);
661 }
662 } else if (nCompReqd == 9) {
663 if (shape==1) {
664 sprintf(tmpBuffer, TENSOR_FORMAT,
665 sampleAvg[0], 0.f, 0.f,
666 0.f, 0.f, 0.f,
667 0.f, 0.f, 0.f);
668 } else if (shape==2) {
669 sprintf(tmpBuffer, TENSOR_FORMAT,
670 sampleAvg[0], sampleAvg[1], 0.f,
671 sampleAvg[2], sampleAvg[3], 0.f,
672 0.f, 0.f, 0.f);
673 } else if (shape==3) {
674 sprintf(tmpBuffer, TENSOR_FORMAT,
675 sampleAvg[0],sampleAvg[1],sampleAvg[2],
676 sampleAvg[3],sampleAvg[4],sampleAvg[5],
677 sampleAvg[6],sampleAvg[7],sampleAvg[8]);
678 }
679 }
680 if ( mpi_size > 1) {
681 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
682 } else {
683 fputs(tmpBuffer, fileHandle_p);
684 }
685 } /* for l (numCellFactor) */
686 } /* if I am the owner */
687 } /* for i (numCells) */
688
689 if ( mpi_size > 1) {
690 MPI_WRITE_ORDERED(txtBuffer);
691 MPI_RANK0_WRITE_SHARED(tag_End_DataArray);
692 } else {
693 fputs(tag_End_DataArray, fileHandle_p);
694 }
695 } /* !isEmpty && cellCentered */
696 } /* for dataIdx */
697
698 strcpy(txtBuffer, "</CellData>\n");
699 if ( mpi_size > 1) {
700 MPI_RANK0_WRITE_SHARED(txtBuffer);
701 } else {
702 fputs(txtBuffer, fileHandle_p);
703 }
704 } /* if noError && writeCellData */
705
706 /************************************************************************/
707 /* write point data */
708
709 if (writePointData && Dudley_noError()) {
710 /* mark the active data arrays */
711 bool_t set_scalar=FALSE, set_vector=FALSE, set_tensor=FALSE;
712 strcpy(txtBuffer, "<PointData");
713 for (dataIdx = 0; dataIdx<num_data; dataIdx++) {
714 if (!isEmpty(data_pp[dataIdx]) && !isCellCentered[dataIdx]) {
715 switch (getDataPointRank(data_pp[dataIdx])) {
716 case 0:
717 if (!set_scalar) {
718 strcat(txtBuffer, " Scalars=\"");
719 strcat(txtBuffer, names_p[dataIdx]);
720 strcat(txtBuffer, "\"");
721 set_scalar = TRUE;
722 }
723 break;
724 case 1:
725 if (!set_vector) {
726 strcat(txtBuffer, " Vectors=\"");
727 strcat(txtBuffer, names_p[dataIdx]);
728 strcat(txtBuffer, "\"");
729 set_vector = TRUE;
730 }
731 break;
732 case 2:
733 if (!set_tensor) {
734 strcat(txtBuffer, " Tensors=\"");
735 strcat(txtBuffer, names_p[dataIdx]);
736 strcat(txtBuffer, "\"");
737 set_tensor = TRUE;
738 }
739 break;
740 default:
741 sprintf(errorMsg, "saveVTK: data %s: VTK supports data with rank <= 2 only.", names_p[dataIdx]);
742 Dudley_setError(VALUE_ERROR, errorMsg);
743 }
744 }
745 if (!Dudley_noError())
746 break;
747 }
748 }
749 /* only continue if no error occurred */
750 if (writePointData && Dudley_noError()) {
751 strcat(txtBuffer, ">\n");
752 if ( mpi_size > 1) {
753 MPI_RANK0_WRITE_SHARED(txtBuffer);
754 } else {
755 fputs(txtBuffer, fileHandle_p);
756 }
757
758 /* write the arrays */
759 for (dataIdx = 0; dataIdx < num_data; dataIdx++) {
760 if (!isEmpty(data_pp[dataIdx]) && !isCellCentered[dataIdx]) {
761 Dudley_NodeMapping* nodeMapping;
762 dim_t rank = getDataPointRank(data_pp[dataIdx]);
763 dim_t nCompReqd = 1; /* number of components mpi_required */
764 if (getFunctionSpaceType(data_pp[dataIdx]) == DUDLEY_REDUCED_NODES) {
765 nodeMapping = mesh_p->Nodes->reducedNodesMapping;
766 } else {
767 nodeMapping = mesh_p->Nodes->nodesMapping;
768 }
769 if (rank == 0) {
770 nCompReqd = 1;
771 shape = 0;
772 } else if (rank == 1) {
773 nCompReqd = 3;
774 shape = getDataPointShape(data_pp[dataIdx], 0);
775 if (shape > 3) {
776 Dudley_setError(VALUE_ERROR, "saveVTK: rank 1 objects must have 3 components at most.");
777 }
778 } else {
779 nCompReqd = 9;
780 shape=getDataPointShape(data_pp[dataIdx], 0);
781 if (shape > 3 || shape != getDataPointShape(data_pp[dataIdx], 1)) {
782 Dudley_setError(VALUE_ERROR, "saveVTK: rank 2 objects of shape 2x2 or 3x3 supported only.");
783 }
784 }
785 /* bail out if an error occurred */
786 if (!Dudley_noError())
787 break;
788
789 sprintf(txtBuffer, tag_Float_DataArray, names_p[dataIdx], nCompReqd);
790 if ( mpi_size > 1) {
791 MPI_RANK0_WRITE_SHARED(txtBuffer);
792 } else {
793 fputs(txtBuffer, fileHandle_p);
794 }
795
796 txtBuffer[0] = '\0';
797 txtBufferInUse = 0;
798 for (i=0; i<mesh_p->Nodes->numNodes; i++) {
799 k = globalNodeIndex[i];
800 if ( (myFirstNode <= k) && (k < myLastNode) ) {
801 __const double *values = getSampleDataRO(data_pp[dataIdx], nodeMapping->target[i]);
802 /* if the number of mpi_required components is more than
803 * the number of actual components, pad with zeros.
804 * Probably only need to get shape of first element */
805 if (nCompReqd == 1) {
806 sprintf(tmpBuffer, SCALAR_FORMAT, values[0]);
807 } else if (nCompReqd == 3) {
808 if (shape==1) {
809 sprintf(tmpBuffer, VECTOR_FORMAT,
810 values[0], 0.f, 0.f);
811 } else if (shape==2) {
812 sprintf(tmpBuffer, VECTOR_FORMAT,
813 values[0], values[1], 0.f);
814 } else if (shape==3) {
815 sprintf(tmpBuffer, VECTOR_FORMAT,
816 values[0], values[1], values[2]);
817 }
818 } else if (nCompReqd == 9) {
819 if (shape==1) {
820 sprintf(tmpBuffer, TENSOR_FORMAT,
821 values[0], 0.f, 0.f,
822 0.f, 0.f, 0.f,
823 0.f, 0.f, 0.f);
824 } else if (shape==2) {
825 sprintf(tmpBuffer, TENSOR_FORMAT,
826 values[0], values[1], 0.f,
827 values[2], values[3], 0.f,
828 0.f, 0.f, 0.f);
829 } else if (shape==3) {
830 sprintf(tmpBuffer, TENSOR_FORMAT,
831 values[0], values[1], values[2],
832 values[3], values[4], values[5],
833 values[6], values[7], values[8]);
834 }
835 }
836 if ( mpi_size > 1) {
837 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
838 } else {
839 fputs(tmpBuffer, fileHandle_p);
840 }
841 } /* if this is my node */
842 } /* for i (numNodes) */
843
844 if ( mpi_size > 1) {
845 MPI_WRITE_ORDERED(txtBuffer);
846 MPI_RANK0_WRITE_SHARED(tag_End_DataArray);
847 } else {
848 fputs(tag_End_DataArray, fileHandle_p);
849 }
850 } /* !isEmpty && !isCellCentered */
851 } /* for dataIdx */
852
853 strcpy(txtBuffer, "</PointData>\n");
854 if ( mpi_size > 1) {
855 MPI_RANK0_WRITE_SHARED(txtBuffer);
856 } else {
857 fputs(txtBuffer, fileHandle_p);
858 }
859 } /* if noError && writePointData */
860
861 /* Final write to VTK file */
862 if (Dudley_noError()) {
863 if (mpi_size > 1) {
864 MPI_RANK0_WRITE_SHARED(vtkFooter);
865 } else {
866 fputs(vtkFooter, fileHandle_p);
867 }
868 }
869
870 if ( mpi_size > 1) {
871 #ifdef PASO_MPI
872 MPI_File_close(&mpi_fileHandle_p);
873 MPI_Barrier(mesh_p->Nodes->MPIInfo->comm);
874 #endif
875 } else {
876 fclose(fileHandle_p);
877 }
878 TMPMEMFREE(isCellCentered);
879 TMPMEMFREE(txtBuffer);
880 }
881

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26