/[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 3086 - (show annotations)
Thu Aug 5 05:07:58 2010 UTC (8 years, 9 months ago) by jfenwick
File MIME type: text/plain
File size: 46445 byte(s)
Another pass at removing finley

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26