/[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 3114 - (show annotations)
Fri Aug 27 05:26:25 2010 UTC (9 years, 6 months ago) by jfenwick
File MIME type: text/plain
File size: 45330 byte(s)
It doesn't pass all tests but this is major progress

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
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 Dudley_Mesh_saveVTK(const char *filename_p,
150 Dudley_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=DUDLEY_NODES, elementType=DUDLEY_UNKNOWN;
179 Dudley_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 Dudley_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 Dudley_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 Dudley_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 Dudley_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 (!Dudley_checkPtr(isCellCentered)) {
264 for (dataIdx=0; dataIdx<num_data; ++dataIdx) {
265 if (! isEmpty(data_pp[dataIdx])) {
266 switch(getFunctionSpaceType(data_pp[dataIdx]) ) {
267 case DUDLEY_NODES:
268 nodeType = (nodeType == DUDLEY_REDUCED_NODES) ? DUDLEY_REDUCED_NODES : DUDLEY_NODES;
269 isCellCentered[dataIdx] = FALSE;
270 if (elementType==DUDLEY_UNKNOWN || elementType==DUDLEY_ELEMENTS) {
271 elementType = DUDLEY_ELEMENTS;
272 } else {
273 Dudley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
274 }
275 break;
276 case DUDLEY_REDUCED_NODES:
277 nodeType = DUDLEY_REDUCED_NODES;
278 isCellCentered[dataIdx] = FALSE;
279 if (elementType==DUDLEY_UNKNOWN || elementType==DUDLEY_ELEMENTS) {
280 elementType = DUDLEY_ELEMENTS;
281 } else {
282 Dudley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
283 }
284 break;
285 case DUDLEY_REDUCED_ELEMENTS:
286 hasReducedElements = TRUE;
287 case DUDLEY_ELEMENTS:
288 isCellCentered[dataIdx] = TRUE;
289 if (elementType==DUDLEY_UNKNOWN || elementType==DUDLEY_ELEMENTS) {
290 elementType = DUDLEY_ELEMENTS;
291 } else {
292 Dudley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
293 }
294 break;
295 case DUDLEY_REDUCED_FACE_ELEMENTS:
296 hasReducedElements = TRUE;
297 case DUDLEY_FACE_ELEMENTS:
298 isCellCentered[dataIdx] = TRUE;
299 if (elementType==DUDLEY_UNKNOWN || elementType==DUDLEY_FACE_ELEMENTS) {
300 elementType = DUDLEY_FACE_ELEMENTS;
301 } else {
302 Dudley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
303 }
304 break;
305 case DUDLEY_POINTS:
306 isCellCentered[dataIdx]=TRUE;
307 if (elementType==DUDLEY_UNKNOWN || elementType==DUDLEY_POINTS) {
308 elementType = DUDLEY_POINTS;
309 } else {
310 Dudley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
311 }
312 break;
313 default:
314 sprintf(errorMsg, "saveVTK: unknown function space type %d",getFunctionSpaceType(data_pp[dataIdx]));
315 Dudley_setError(TYPE_ERROR, errorMsg);
316 }
317 if (isCellCentered[dataIdx]) {
318 writeCellData = TRUE;
319 } else {
320 writePointData = TRUE;
321 }
322 maxNameLen = MAX(maxNameLen, strlen(names_p[dataIdx]));
323 }
324 }
325 }
326
327 /************************************************************************/
328 /* select number of points and the mesh component */
329
330 if (Dudley_noError()) {
331 if (nodeType == DUDLEY_REDUCED_NODES) {
332 myFirstNode = Dudley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
333 myLastNode = Dudley_NodeFile_getLastReducedNode(mesh_p->Nodes);
334 globalNumPoints = Dudley_NodeFile_getGlobalNumReducedNodes(mesh_p->Nodes);
335 globalNodeIndex = Dudley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
336 } else {
337 myFirstNode = Dudley_NodeFile_getFirstNode(mesh_p->Nodes);
338 myLastNode = Dudley_NodeFile_getLastNode(mesh_p->Nodes);
339 globalNumPoints = Dudley_NodeFile_getGlobalNumNodes(mesh_p->Nodes);
340 globalNodeIndex = Dudley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
341 }
342 myNumPoints = myLastNode - myFirstNode;
343 if (elementType==DUDLEY_UNKNOWN) elementType=DUDLEY_ELEMENTS;
344 switch(elementType) {
345 case DUDLEY_ELEMENTS:
346 elements = mesh_p->Elements;
347 break;
348 case DUDLEY_FACE_ELEMENTS:
349 elements = mesh_p->FaceElements;
350 break;
351 case DUDLEY_POINTS:
352 elements = mesh_p->Points;
353 break;
354 }
355 if (elements==NULL) {
356 Dudley_setError(SYSTEM_ERROR, "saveVTK: undefined element file");
357 } else {
358 /* map dudley element type to VTK element type */
359 numCells = elements->numElements;
360 globalNumCells = Dudley_ElementFile_getGlobalNumElements(elements);
361 myNumCells = Dudley_ElementFile_getMyNumElements(elements);
362 myFirstCell = Dudley_ElementFile_getFirstElement(elements);
363 NN = elements->numNodes;
364 if (hasReducedElements) {
365 quadNodes_p=elements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->QuadNodes;
366 } else {
367 quadNodes_p=elements->referenceElementSet->referenceElement->Parametrization->QuadNodes;
368 }
369 if (nodeType==DUDLEY_REDUCED_NODES) {
370 typeId = elements->referenceElementSet->referenceElement->LinearType->TypeId;
371 } else {
372 typeId = elements->referenceElementSet->referenceElement->Type->TypeId;
373 }
374 switch (typeId) {
375 case Point1:
376 case Line2Face:
377 case Line3Face:
378 case Point1_Contact:
379 case Line2Face_Contact:
380 case Line3Face_Contact:
381 cellType = VTK_VERTEX;
382 numVTKNodesPerElement = 1;
383 break;
384
385 case Line2:
386 case Tri3Face:
387 // case Rec4Face:
388 case Line2_Contact:
389 case Tri3_Contact:
390 case Tri3Face_Contact:
391 /* case Rec4Face_Contact:*/
392 cellType = VTK_LINE;
393 numVTKNodesPerElement = 2;
394 break;
395
396 case Line3Macro:
397 cellType = VTK_LINE;
398 numCellFactor = 2;
399 numVTKNodesPerElement = 2;
400 break;
401
402 case Tri3:
403 case Tet4Face:
404 case Tet4Face_Contact:
405 cellType = VTK_TRIANGLE;
406 numVTKNodesPerElement = 3;
407 break;
408
409 // case Rec4:
410 // case Hex8Face:
411 // case Rec4_Contact:
412 // case Hex8Face_Contact:
413 // cellType = VTK_QUAD;
414 // numVTKNodesPerElement = 4;
415 // break;
416
417 // case Rec9Macro:
418 // case Rec9:
419 // numCellFactor = 4;
420 // cellType = VTK_QUAD;
421 // numVTKNodesPerElement = 4;
422 // break;
423
424 case Tet4:
425 cellType = VTK_TETRA;
426 numVTKNodesPerElement = 4;
427 break;
428
429 // case Hex8:
430 // cellType = VTK_HEXAHEDRON;
431 // numVTKNodesPerElement = 8;
432 // break;
433
434 case Line3:
435 case Tri6Face:
436 // case Rec8Face:
437 case Line3_Contact:
438 case Tri6Face_Contact:
439 // case Rec8Face_Contact:
440 cellType = VTK_QUADRATIC_EDGE;
441 numVTKNodesPerElement = 3;
442 break;
443
444 case Tri6:
445 case Tri6Macro:
446 case Tet10Face:
447 case Tri6_Contact:
448 case Tet10Face_Contact:
449 cellType = VTK_QUADRATIC_TRIANGLE;
450 numVTKNodesPerElement = 6;
451 break;
452
453
454 // case Rec8:
455 // case Hex20Face:
456 // case Rec8_Contact:
457 // case Hex20Face_Contact:
458 // cellType = VTK_QUADRATIC_QUAD;
459 // numVTKNodesPerElement = 8;
460 // break;
461
462 case Tet10Macro:
463 case Tet10:
464 cellType = VTK_QUADRATIC_TETRA;
465 numVTKNodesPerElement = 10;
466 break;
467
468 // case Hex20:
469 // cellType = VTK_QUADRATIC_HEXAHEDRON;
470 // numVTKNodesPerElement = 20;
471 // break;
472
473 // case Hex27Macro:
474 // case Hex27:
475 // numCellFactor = 8;
476 // cellType = VTK_HEXAHEDRON;
477 // numVTKNodesPerElement = 8;
478 // break;
479
480 default:
481 sprintf(errorMsg, "saveVTK: Element type %s is not supported by VTK.", elements->referenceElementSet->referenceElement->Type->Name);
482 Dudley_setError(VALUE_ERROR, errorMsg);
483 }
484 }
485 }
486
487 /* allocate enough memory for text buffer */
488
489 txtBufferSize = strlen(vtkHeader) + 3*LEN_INT_FORMAT + (30+3*maxNameLen)+strlen(metadata)+strlen(metadata_schema);
490 if (mpi_size > 1) {
491 txtBufferSize = MAX(txtBufferSize, myNumPoints * LEN_TMP_BUFFER);
492 txtBufferSize = MAX(txtBufferSize, numCellFactor * myNumCells *
493 (LEN_INT_FORMAT * numVTKNodesPerElement + 1));
494 txtBufferSize = MAX(txtBufferSize,
495 numCellFactor * myNumCells * LEN_TENSOR_FORMAT);
496 txtBufferSize = MAX(txtBufferSize, myNumPoints * LEN_TENSOR_FORMAT);
497 }
498 txtBuffer = TMPMEMALLOC(txtBufferSize+1, char);
499
500 /* sets error if memory allocation failed */
501 Dudley_checkPtr(txtBuffer);
502
503 /************************************************************************/
504 /* write number of points and the mesh component */
505
506 if (Dudley_noError()) {
507 const index_t *nodeIndex;
508 if (DUDLEY_REDUCED_NODES == nodeType) {
509 nodeIndex = elements->referenceElementSet->referenceElement->Type->linearNodes;
510 } else if (Line3Macro == typeId) {
511 nodeIndex = VTK_LINE3_INDEX;
512 }/* else if ( (Rec9 == typeId) || (Rec9Macro == typeId) ) {
513 nodeIndex = VTK_REC9_INDEX;
514 } else if (Hex20 == typeId) {
515 nodeIndex = VTK_HEX20_INDEX;
516 } else if ( (Hex27 == typeId) || (Hex27Macro == typeId) ){
517 nodeIndex = VTK_HEX27_INDEX;
518 } */ else if (numVTKNodesPerElement != elements->referenceElementSet->referenceElement->Type->numNodes) {
519 nodeIndex = elements->referenceElementSet->referenceElement->Type->relevantGeoNodes;
520 } else {
521 nodeIndex = NULL;
522 }
523
524 if (strlen(metadata)>0) {
525 if (strlen(metadata_schema)>0) {
526 sprintf(txtBuffer, vtkHeader," ",metadata_schema,metadata,"\n",globalNumPoints, numCellFactor*globalNumCells, 3);
527 } else {
528 sprintf(txtBuffer, vtkHeader,"","",metadata,"\n",globalNumPoints, numCellFactor*globalNumCells, 3);
529 }
530 } else {
531 if (strlen(metadata_schema)>0) {
532 sprintf(txtBuffer, vtkHeader," ",metadata_schema,"","",globalNumPoints, numCellFactor*globalNumCells, 3);
533 } else {
534 sprintf(txtBuffer, vtkHeader,"","","","",globalNumPoints, numCellFactor*globalNumCells, 3);
535 }
536 }
537
538 if (mpi_size > 1) {
539 /* write the nodes */
540 MPI_RANK0_WRITE_SHARED(txtBuffer);
541 txtBuffer[0] = '\0';
542 txtBufferInUse = 0;
543 if (nDim==2) {
544 for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
545 if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
546 sprintf(tmpBuffer, VECTOR_FORMAT,
547 mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
548 mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
549 0.);
550 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
551 }
552 }
553 } else {
554 for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
555 if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
556 sprintf(tmpBuffer, VECTOR_FORMAT,
557 mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
558 mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
559 mesh_p->Nodes->Coordinates[INDEX2(2, i, nDim)]);
560 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
561 }
562 }
563 } /* nDim */
564 MPI_WRITE_ORDERED(txtBuffer);
565
566 /* write the cells */
567 MPI_RANK0_WRITE_SHARED(tags_End_Points_and_Start_Conn);
568 txtBuffer[0] = '\0';
569 txtBufferInUse = 0;
570 if (nodeIndex == NULL) {
571 for (i = 0; i < numCells; i++) {
572 if (elements->Owner[i] == my_mpi_rank) {
573 for (j = 0; j < numVTKNodesPerElement; j++) {
574 sprintf(tmpBuffer, INT_FORMAT, globalNodeIndex[elements->Nodes[INDEX2(j, i, NN)]]);
575 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
576 }
577 __STRCAT(txtBuffer, NEWLINE, txtBufferInUse);
578 }
579 }
580 } else {
581 for (i = 0; i < numCells; i++) {
582 if (elements->Owner[i] == my_mpi_rank) {
583 for (l = 0; l < numCellFactor; l++) {
584 const int* idx=&nodeIndex[l*numVTKNodesPerElement];
585 for (j = 0; j < numVTKNodesPerElement; j++) {
586 sprintf(tmpBuffer, INT_FORMAT, globalNodeIndex[elements->Nodes[INDEX2(idx[j], i, NN)]]);
587 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
588 }
589 __STRCAT(txtBuffer, NEWLINE, txtBufferInUse);
590 }
591 }
592 }
593 } /* nodeIndex */
594 MPI_WRITE_ORDERED(txtBuffer);
595
596 /* write the offsets */
597 MPI_RANK0_WRITE_SHARED(tags_End_Conn_and_Start_Offset);
598 txtBuffer[0] = '\0';
599 txtBufferInUse = 0;
600 for (i = numVTKNodesPerElement*(myFirstCell*numCellFactor+1);
601 i <= (myFirstCell+myNumCells)*numVTKNodesPerElement*numCellFactor; i += numVTKNodesPerElement)
602 {
603 sprintf(tmpBuffer, INT_NEWLINE_FORMAT, i);
604 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
605 }
606 MPI_WRITE_ORDERED(txtBuffer);
607
608 /* write element type */
609 sprintf(tmpBuffer, INT_NEWLINE_FORMAT, cellType);
610 MPI_RANK0_WRITE_SHARED(tags_End_Offset_and_Start_Type);
611 txtBuffer[0] = '\0';
612 txtBufferInUse = 0;
613 for (i = numVTKNodesPerElement*(myFirstCell*numCellFactor+1);
614 i <= (myFirstCell+myNumCells)*numVTKNodesPerElement*numCellFactor; i += numVTKNodesPerElement)
615 {
616 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
617 }
618 MPI_WRITE_ORDERED(txtBuffer);
619 /* finalize cell information */
620 strcpy(txtBuffer, "</DataArray>\n</Cells>\n");
621 MPI_RANK0_WRITE_SHARED(txtBuffer);
622
623 } else { /***** mpi_size == 1 *****/
624
625 /* write the nodes */
626 fputs(txtBuffer, fileHandle_p);
627 if (nDim==2) {
628 for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
629 if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
630 fprintf(fileHandle_p, VECTOR_FORMAT,
631 mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
632 mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
633 0.);
634 }
635 }
636 } else {
637 for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
638 if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
639 fprintf(fileHandle_p, VECTOR_FORMAT,
640 mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
641 mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
642 mesh_p->Nodes->Coordinates[INDEX2(2, i, nDim)]);
643 }
644 }
645 } /* nDim */
646
647 /* write the cells */
648 fputs(tags_End_Points_and_Start_Conn, fileHandle_p);
649 if (nodeIndex == NULL) {
650 for (i = 0; i < numCells; i++) {
651 for (j = 0; j < numVTKNodesPerElement; j++) {
652 fprintf(fileHandle_p, INT_FORMAT, globalNodeIndex[elements->Nodes[INDEX2(j, i, NN)]]);
653 }
654 fprintf(fileHandle_p, NEWLINE);
655 }
656 } else {
657 for (i = 0; i < numCells; i++) {
658 for (l = 0; l < numCellFactor; l++) {
659 const int* idx=&nodeIndex[l*numVTKNodesPerElement];
660 for (j = 0; j < numVTKNodesPerElement; j++) {
661 fprintf(fileHandle_p, INT_FORMAT, globalNodeIndex[elements->Nodes[INDEX2(idx[j], i, NN)]]);
662 }
663 fprintf(fileHandle_p, NEWLINE);
664 }
665 }
666 } /* nodeIndex */
667
668 /* write the offsets */
669 fputs(tags_End_Conn_and_Start_Offset, fileHandle_p);
670 for (i = numVTKNodesPerElement; i <= numCells*numVTKNodesPerElement*numCellFactor; i += numVTKNodesPerElement) {
671 fprintf(fileHandle_p, INT_NEWLINE_FORMAT, i);
672 }
673
674 /* write element type */
675 sprintf(tmpBuffer, INT_NEWLINE_FORMAT, cellType);
676 fputs(tags_End_Offset_and_Start_Type, fileHandle_p);
677 for (i = 0; i < numCells*numCellFactor; i++)
678 fputs(tmpBuffer, fileHandle_p);
679 /* finalize cell information */
680 fputs("</DataArray>\n</Cells>\n", fileHandle_p);
681 } /* mpi_size */
682
683 } /* Dudley_noError */
684
685 /************************************************************************/
686 /* write cell data */
687
688 if (writeCellData && Dudley_noError()) {
689 bool_t set_scalar=FALSE, set_vector=FALSE, set_tensor=FALSE;
690 /* mark the active data arrays */
691 strcpy(txtBuffer, "<CellData");
692 for (dataIdx = 0; dataIdx < num_data; dataIdx++) {
693 if (!isEmpty(data_pp[dataIdx]) && isCellCentered[dataIdx]) {
694 /* rank == 0 <--> scalar data */
695 /* rank == 1 <--> vector data */
696 /* rank == 2 <--> tensor data */
697 switch (getDataPointRank(data_pp[dataIdx])) {
698 case 0:
699 if (!set_scalar) {
700 strcat(txtBuffer, " Scalars=\"");
701 strcat(txtBuffer, names_p[dataIdx]);
702 strcat(txtBuffer, "\"");
703 set_scalar = TRUE;
704 }
705 break;
706 case 1:
707 if (!set_vector) {
708 strcat(txtBuffer, " Vectors=\"");
709 strcat(txtBuffer, names_p[dataIdx]);
710 strcat(txtBuffer, "\"");
711 set_vector = TRUE;
712 }
713 break;
714 case 2:
715 if (!set_tensor) {
716 strcat(txtBuffer, " Tensors=\"");
717 strcat(txtBuffer, names_p[dataIdx]);
718 strcat(txtBuffer, "\"");
719 set_tensor = TRUE;
720 }
721 break;
722 default:
723 sprintf(errorMsg, "saveVTK: data %s: VTK supports data with rank <= 2 only.", names_p[dataIdx]);
724 Dudley_setError(VALUE_ERROR, errorMsg);
725 }
726 }
727 if (!Dudley_noError())
728 break;
729 }
730 }
731 /* only continue if no error occurred */
732 if (writeCellData && Dudley_noError()) {
733 strcat(txtBuffer, ">\n");
734 if ( mpi_size > 1) {
735 MPI_RANK0_WRITE_SHARED(txtBuffer);
736 } else {
737 fputs(txtBuffer, fileHandle_p);
738 }
739
740 /* write the arrays */
741 for (dataIdx = 0; dataIdx < num_data; dataIdx++) {
742 if (!isEmpty(data_pp[dataIdx]) && isCellCentered[dataIdx]) {
743 dim_t numPointsPerSample=getNumDataPointsPerSample(data_pp[dataIdx]);
744 dim_t rank = getDataPointRank(data_pp[dataIdx]);
745 dim_t nComp = getDataPointSize(data_pp[dataIdx]);
746 dim_t nCompReqd = 1; /* number of components mpi_required */
747 if (rank == 0) {
748 nCompReqd = 1;
749 shape = 0;
750 } else if (rank == 1) {
751 nCompReqd = 3;
752 shape = getDataPointShape(data_pp[dataIdx], 0);
753 if (shape > 3) {
754 Dudley_setError(VALUE_ERROR, "saveVTK: rank 1 objects must have 3 components at most.");
755 }
756 } else {
757 nCompReqd = 9;
758 shape = getDataPointShape(data_pp[dataIdx], 0);
759 if (shape > 3 || shape != getDataPointShape(data_pp[dataIdx], 1)) {
760 Dudley_setError(VALUE_ERROR, "saveVTK: rank 2 objects of shape 2x2 or 3x3 supported only.");
761 }
762 }
763 /* bail out if an error occurred */
764 if (!Dudley_noError())
765 break;
766
767 sprintf(txtBuffer, tag_Float_DataArray, names_p[dataIdx], nCompReqd);
768 if ( mpi_size > 1) {
769 MPI_RANK0_WRITE_SHARED(txtBuffer);
770 } else {
771 fputs(txtBuffer, fileHandle_p);
772 }
773
774 txtBuffer[0] = '\0';
775 txtBufferInUse = 0;
776 for (i=0; i<numCells; i++) {
777 if (elements->Owner[i] == my_mpi_rank) {
778 __const double *values = getSampleDataRO(data_pp[dataIdx], i);
779 for (l = 0; l < numCellFactor; l++) {
780 double sampleAvg[NCOMP_MAX];
781 dim_t nCompUsed = MIN(nComp, NCOMP_MAX);
782
783 /* average over number of points in the sample */
784 if (isExpanded(data_pp[dataIdx])) {
785 dim_t hits=0;
786 for (k=0; k<nCompUsed; k++) sampleAvg[k]=0;
787 for (j=0; j<numPointsPerSample; j++) {
788 if (nodeInQuadrant(quadNodes_p, typeId, j, l)) {
789 hits++;
790 for (k=0; k<nCompUsed; k++) {
791 sampleAvg[k] += values[INDEX2(k,j,nComp)];
792 }
793 }
794 }
795 for (k=0; k<nCompUsed; k++)
796 sampleAvg[k] /= MAX(hits, 1);
797 } else {
798 for (k=0; k<nCompUsed; k++)
799 sampleAvg[k] = values[k];
800 } /* isExpanded */
801
802 /* if the number of required components is more than
803 * the number of actual components, pad with zeros
804 */
805 /* probably only need to get shape of first element */
806 if (nCompReqd == 1) {
807 sprintf(tmpBuffer, SCALAR_FORMAT, sampleAvg[0]);
808 } else if (nCompReqd == 3) {
809 if (shape==1) {
810 sprintf(tmpBuffer, VECTOR_FORMAT,
811 sampleAvg[0], 0.f, 0.f);
812 } else if (shape==2) {
813 sprintf(tmpBuffer, VECTOR_FORMAT,
814 sampleAvg[0], sampleAvg[1], 0.f);
815 } else if (shape==3) {
816 sprintf(tmpBuffer, VECTOR_FORMAT,
817 sampleAvg[0],sampleAvg[1],sampleAvg[2]);
818 }
819 } else if (nCompReqd == 9) {
820 if (shape==1) {
821 sprintf(tmpBuffer, TENSOR_FORMAT,
822 sampleAvg[0], 0.f, 0.f,
823 0.f, 0.f, 0.f,
824 0.f, 0.f, 0.f);
825 } else if (shape==2) {
826 sprintf(tmpBuffer, TENSOR_FORMAT,
827 sampleAvg[0], sampleAvg[1], 0.f,
828 sampleAvg[2], sampleAvg[3], 0.f,
829 0.f, 0.f, 0.f);
830 } else if (shape==3) {
831 sprintf(tmpBuffer, TENSOR_FORMAT,
832 sampleAvg[0],sampleAvg[1],sampleAvg[2],
833 sampleAvg[3],sampleAvg[4],sampleAvg[5],
834 sampleAvg[6],sampleAvg[7],sampleAvg[8]);
835 }
836 }
837 if ( mpi_size > 1) {
838 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
839 } else {
840 fputs(tmpBuffer, fileHandle_p);
841 }
842 } /* for l (numCellFactor) */
843 } /* if I am the owner */
844 } /* for i (numCells) */
845
846 if ( mpi_size > 1) {
847 MPI_WRITE_ORDERED(txtBuffer);
848 MPI_RANK0_WRITE_SHARED(tag_End_DataArray);
849 } else {
850 fputs(tag_End_DataArray, fileHandle_p);
851 }
852 } /* !isEmpty && cellCentered */
853 } /* for dataIdx */
854
855 strcpy(txtBuffer, "</CellData>\n");
856 if ( mpi_size > 1) {
857 MPI_RANK0_WRITE_SHARED(txtBuffer);
858 } else {
859 fputs(txtBuffer, fileHandle_p);
860 }
861 } /* if noError && writeCellData */
862
863 /************************************************************************/
864 /* write point data */
865
866 if (writePointData && Dudley_noError()) {
867 /* mark the active data arrays */
868 bool_t set_scalar=FALSE, set_vector=FALSE, set_tensor=FALSE;
869 strcpy(txtBuffer, "<PointData");
870 for (dataIdx = 0; dataIdx<num_data; dataIdx++) {
871 if (!isEmpty(data_pp[dataIdx]) && !isCellCentered[dataIdx]) {
872 switch (getDataPointRank(data_pp[dataIdx])) {
873 case 0:
874 if (!set_scalar) {
875 strcat(txtBuffer, " Scalars=\"");
876 strcat(txtBuffer, names_p[dataIdx]);
877 strcat(txtBuffer, "\"");
878 set_scalar = TRUE;
879 }
880 break;
881 case 1:
882 if (!set_vector) {
883 strcat(txtBuffer, " Vectors=\"");
884 strcat(txtBuffer, names_p[dataIdx]);
885 strcat(txtBuffer, "\"");
886 set_vector = TRUE;
887 }
888 break;
889 case 2:
890 if (!set_tensor) {
891 strcat(txtBuffer, " Tensors=\"");
892 strcat(txtBuffer, names_p[dataIdx]);
893 strcat(txtBuffer, "\"");
894 set_tensor = TRUE;
895 }
896 break;
897 default:
898 sprintf(errorMsg, "saveVTK: data %s: VTK supports data with rank <= 2 only.", names_p[dataIdx]);
899 Dudley_setError(VALUE_ERROR, errorMsg);
900 }
901 }
902 if (!Dudley_noError())
903 break;
904 }
905 }
906 /* only continue if no error occurred */
907 if (writePointData && Dudley_noError()) {
908 strcat(txtBuffer, ">\n");
909 if ( mpi_size > 1) {
910 MPI_RANK0_WRITE_SHARED(txtBuffer);
911 } else {
912 fputs(txtBuffer, fileHandle_p);
913 }
914
915 /* write the arrays */
916 for (dataIdx = 0; dataIdx < num_data; dataIdx++) {
917 if (!isEmpty(data_pp[dataIdx]) && !isCellCentered[dataIdx]) {
918 Dudley_NodeMapping* nodeMapping;
919 dim_t rank = getDataPointRank(data_pp[dataIdx]);
920 dim_t nCompReqd = 1; /* number of components mpi_required */
921 if (getFunctionSpaceType(data_pp[dataIdx]) == DUDLEY_REDUCED_NODES) {
922 nodeMapping = mesh_p->Nodes->reducedNodesMapping;
923 } else {
924 nodeMapping = mesh_p->Nodes->nodesMapping;
925 }
926 if (rank == 0) {
927 nCompReqd = 1;
928 shape = 0;
929 } else if (rank == 1) {
930 nCompReqd = 3;
931 shape = getDataPointShape(data_pp[dataIdx], 0);
932 if (shape > 3) {
933 Dudley_setError(VALUE_ERROR, "saveVTK: rank 1 objects must have 3 components at most.");
934 }
935 } else {
936 nCompReqd = 9;
937 shape=getDataPointShape(data_pp[dataIdx], 0);
938 if (shape > 3 || shape != getDataPointShape(data_pp[dataIdx], 1)) {
939 Dudley_setError(VALUE_ERROR, "saveVTK: rank 2 objects of shape 2x2 or 3x3 supported only.");
940 }
941 }
942 /* bail out if an error occurred */
943 if (!Dudley_noError())
944 break;
945
946 sprintf(txtBuffer, tag_Float_DataArray, names_p[dataIdx], nCompReqd);
947 if ( mpi_size > 1) {
948 MPI_RANK0_WRITE_SHARED(txtBuffer);
949 } else {
950 fputs(txtBuffer, fileHandle_p);
951 }
952
953 txtBuffer[0] = '\0';
954 txtBufferInUse = 0;
955 for (i=0; i<mesh_p->Nodes->numNodes; i++) {
956 k = globalNodeIndex[i];
957 if ( (myFirstNode <= k) && (k < myLastNode) ) {
958 __const double *values = getSampleDataRO(data_pp[dataIdx], nodeMapping->target[i]);
959 /* if the number of mpi_required components is more than
960 * the number of actual components, pad with zeros.
961 * Probably only need to get shape of first element */
962 if (nCompReqd == 1) {
963 sprintf(tmpBuffer, SCALAR_FORMAT, values[0]);
964 } else if (nCompReqd == 3) {
965 if (shape==1) {
966 sprintf(tmpBuffer, VECTOR_FORMAT,
967 values[0], 0.f, 0.f);
968 } else if (shape==2) {
969 sprintf(tmpBuffer, VECTOR_FORMAT,
970 values[0], values[1], 0.f);
971 } else if (shape==3) {
972 sprintf(tmpBuffer, VECTOR_FORMAT,
973 values[0], values[1], values[2]);
974 }
975 } else if (nCompReqd == 9) {
976 if (shape==1) {
977 sprintf(tmpBuffer, TENSOR_FORMAT,
978 values[0], 0.f, 0.f,
979 0.f, 0.f, 0.f,
980 0.f, 0.f, 0.f);
981 } else if (shape==2) {
982 sprintf(tmpBuffer, TENSOR_FORMAT,
983 values[0], values[1], 0.f,
984 values[2], values[3], 0.f,
985 0.f, 0.f, 0.f);
986 } else if (shape==3) {
987 sprintf(tmpBuffer, TENSOR_FORMAT,
988 values[0], values[1], values[2],
989 values[3], values[4], values[5],
990 values[6], values[7], values[8]);
991 }
992 }
993 if ( mpi_size > 1) {
994 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
995 } else {
996 fputs(tmpBuffer, fileHandle_p);
997 }
998 } /* if this is my node */
999 } /* for i (numNodes) */
1000
1001 if ( mpi_size > 1) {
1002 MPI_WRITE_ORDERED(txtBuffer);
1003 MPI_RANK0_WRITE_SHARED(tag_End_DataArray);
1004 } else {
1005 fputs(tag_End_DataArray, fileHandle_p);
1006 }
1007 } /* !isEmpty && !isCellCentered */
1008 } /* for dataIdx */
1009
1010 strcpy(txtBuffer, "</PointData>\n");
1011 if ( mpi_size > 1) {
1012 MPI_RANK0_WRITE_SHARED(txtBuffer);
1013 } else {
1014 fputs(txtBuffer, fileHandle_p);
1015 }
1016 } /* if noError && writePointData */
1017
1018 /* Final write to VTK file */
1019 if (Dudley_noError()) {
1020 if (mpi_size > 1) {
1021 MPI_RANK0_WRITE_SHARED(vtkFooter);
1022 } else {
1023 fputs(vtkFooter, fileHandle_p);
1024 }
1025 }
1026
1027 if ( mpi_size > 1) {
1028 #ifdef PASO_MPI
1029 MPI_File_close(&mpi_fileHandle_p);
1030 MPI_Barrier(mesh_p->Nodes->MPIInfo->comm);
1031 #endif
1032 } else {
1033 fclose(fileHandle_p);
1034 }
1035 TMPMEMFREE(isCellCentered);
1036 TMPMEMFREE(txtBuffer);
1037 }
1038

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26