/[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 3220 - (show annotations)
Wed Sep 29 00:33:16 2010 UTC (8 years, 7 months ago) by jfenwick
File MIME type: text/plain
File size: 39647 byte(s)
Removing references to Quadrature.? and ShapeFns

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26