/[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 2421 - (show annotations)
Thu May 14 02:29:06 2009 UTC (10 years ago) by gross
Original Path: trunk/finley/src/Mesh_saveVTK.c
File MIME type: text/plain
File size: 45983 byte(s)
vtk writer supports now meta data and meta data schema
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2008 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 FINLEY_NODES or FINLEY_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+1)
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 #define LEN_TENSOR_FORMAT (unsigned int)(9*(12+1)+1)
31 #define NEWLINE "\n"
32 #define LEN_TMP_BUFFER LEN_TENSOR_FORMAT+(MAX_numNodes*LEN_INT_FORMAT+1)+2
33 #define NCOMP_MAX (unsigned int)9
34
35 #define __STRCAT(dest, chunk, dest_in_use) \
36 do {\
37 strcpy(&dest[dest_in_use], chunk);\
38 dest_in_use += strlen(chunk);\
39 } while(0)
40
41 #ifdef PASO_MPI
42 /* writes buffer to file catching the empty buffer case which causes problems
43 * with some MPI versions */
44 #define MPI_WRITE_ORDERED(BUF) \
45 do {\
46 int LLEN=0; \
47 LLEN=(int) strlen(BUF); \
48 if (LLEN==0) { strcpy(BUF, ""); LLEN=0; }\
49 MPI_File_write_ordered(mpi_fileHandle_p, BUF, LLEN, MPI_CHAR, &mpi_status);\
50 } while(0)
51
52 /* writes buffer to file on master only */
53 #define MPI_RANK0_WRITE_SHARED(BUF) \
54 do {\
55 int LLEN=0; \
56 if (my_mpi_rank == 0) {\
57 LLEN=(int) strlen(BUF); \
58 if (LLEN==0) { strcpy(BUF,""); LLEN=0; }\
59 MPI_File_iwrite_shared(mpi_fileHandle_p, BUF, LLEN, MPI_CHAR, &mpi_req);\
60 MPI_Wait(&mpi_req, &mpi_status);\
61 }\
62 } while(0)
63
64 /* For reference only. Investigation needed as to which values may improve
65 * performance */
66 #if 0
67 void create_MPIInfo(MPI_Info& info)
68 {
69 MPI_Info_create(&info);
70 MPI_Info_set(info, "access_style", "write_once, sequential");
71 MPI_Info_set(info, "collective_buffering", "true");
72 MPI_Info_set(info, "cb_block_size", "131072");
73 MPI_Info_set(info, "cb_buffer_size", "1048567");
74 MPI_Info_set(info, "cb_nodes", "8");
75 MPI_Info_set(info, "striping_factor", "16");
76 MPI_Info_set(info, "striping_unit", "424288");
77 }
78 #endif
79
80 #else
81
82 #define MPI_WRITE_ORDERED(A)
83 #define MPI_RANK0_WRITE_SHARED(A)
84
85 #endif /* PASO_MPI */
86
87
88 /* Returns one if the node given by coords and idx is within the quadrant
89 * indexed by q and if the element type is Rec9 or Hex27, zero otherwise */
90 int nodeInQuadrant(const double *coords, ElementTypeId type, int idx, int q)
91 {
92 #define INSIDE_1D(_X_,_C_,_R_) ( ABS((_X_)-(_C_)) <= (_R_) )
93 #define INSIDE_2D(_X_,_Y_,_CX_,_CY_,_R_) ( INSIDE_1D(_X_,_CX_,_R_) && INSIDE_1D(_Y_,_CY_,_R_))
94 #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_) )
95
96 int ret;
97 if (type == Rec9) {
98 if (q==0)
99 ret = INSIDE_2D(coords[2*idx], coords[2*idx+1], 0.25, 0.25, 0.25);
100 else if (q==1)
101 ret = INSIDE_2D(coords[2*idx], coords[2*idx+1], 0.75, 0.25, 0.25);
102 else if (q==2)
103 ret = INSIDE_2D(coords[2*idx], coords[2*idx+1], 0.25, 0.75, 0.25);
104 else if (q==3)
105 ret = INSIDE_2D(coords[2*idx], coords[2*idx+1], 0.75, 0.75, 0.25);
106 else
107 ret = 0;
108 } else if (type == Hex27) {
109 if (q==0)
110 ret = INSIDE_3D(coords[3*idx], coords[3*idx+1], coords[3*idx+2],
111 0.25, 0.25, 0.25, 0.25);
112 else if (q==1)
113 ret = INSIDE_3D(coords[3*idx], coords[3*idx+1], coords[3*idx+2],
114 0.75, 0.25, 0.25, 0.25);
115 else if (q==2)
116 ret = INSIDE_3D(coords[3*idx], coords[3*idx+1], coords[3*idx+2],
117 0.25, 0.75, 0.25, 0.25);
118 else if (q==3)
119 ret = INSIDE_3D(coords[3*idx], coords[3*idx+1], coords[3*idx+2],
120 0.75, 0.75, 0.25, 0.25);
121 else if (q==4)
122 ret = INSIDE_3D(coords[3*idx], coords[3*idx+1], coords[3*idx+2],
123 0.25, 0.25, 0.75, 0.25);
124 else if (q==5)
125 ret = INSIDE_3D(coords[3*idx], coords[3*idx+1], coords[3*idx+2],
126 0.75, 0.25, 0.75, 0.25);
127 else if (q==6)
128 ret = INSIDE_3D(coords[3*idx], coords[3*idx+1], coords[3*idx+2],
129 0.25, 0.75, 0.75, 0.25);
130 else if (q==7)
131 ret = INSIDE_3D(coords[3*idx], coords[3*idx+1], coords[3*idx+2],
132 0.75, 0.75, 0.75, 0.25);
133 else
134 ret = 0;
135 } else {
136 ret = 1;
137 }
138 return ret;
139 }
140
141 void Finley_Mesh_saveVTK(const char *filename_p,
142 Finley_Mesh *mesh_p,
143 const dim_t num_data,
144 char **names_p,
145 escriptDataC **data_pp,
146 const char* metadata,
147 const char*metadata_schema)
148 {
149 #ifdef PASO_MPI
150 MPI_File mpi_fileHandle_p;
151 MPI_Status mpi_status;
152 MPI_Request mpi_req;
153 MPI_Info mpi_info = MPI_INFO_NULL;
154 #endif
155 Paso_MPI_rank my_mpi_rank;
156 FILE *fileHandle_p = NULL;
157 char errorMsg[LenErrorMsg_MAX], *txtBuffer;
158 char tmpBuffer[LEN_TMP_BUFFER];
159 size_t txtBufferSize, txtBufferInUse, maxNameLen;
160 double *quadNodes_p = NULL;
161 dim_t dataIdx, nDim;
162 dim_t numCells=0, globalNumCells=0, numVTKNodesPerElement=0;
163 dim_t myNumPoints=0, globalNumPoints=0;
164 dim_t shape, NN=0, numCellFactor=1, myNumCells=0;
165 bool_t *isCellCentered;
166 bool_t writeCellData=FALSE, writePointData=FALSE, hasReducedElements=FALSE;
167 index_t myFirstNode=0, myLastNode=0, *globalNodeIndex=NULL;
168 index_t myFirstCell=0, k;
169 int mpi_size, i, j, l;
170 int cellType=0, nodeType=FINLEY_NODES, elementType=FINLEY_UNKNOWN;
171 Finley_ElementFile *elements = NULL;
172 ElementTypeId typeId = NoType;
173
174 const char *vtkHeader = \
175 "<?xml version=\"1.0\"?>\n" \
176 "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\"%s%s>\n%s%s" \
177 "<UnstructuredGrid>\n" \
178 "<Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n" \
179 "<Points>\n" \
180 "<DataArray NumberOfComponents=\"%d\" type=\"Float64\" format=\"ascii\">\n";
181 char *vtkFooter = "</Piece>\n</UnstructuredGrid>\n</VTKFile>\n";
182 const char *tag_Float_DataArray="<DataArray Name=\"%s\" type=\"Float64\" NumberOfComponents=\"%d\" format=\"ascii\">\n";
183 char *tags_End_Points_and_Start_Conn = "</DataArray>\n</Points>\n<Cells>\n<DataArray Name=\"connectivity\" type=\"Int32\" format=\"ascii\">\n" ;
184 char *tags_End_Conn_and_Start_Offset = "</DataArray>\n<DataArray Name=\"offsets\" type=\"Int32\" format=\"ascii\">\n";
185 char *tags_End_Offset_and_Start_Type = "</DataArray>\n<DataArray Name=\"types\" type=\"UInt8\" format=\"ascii\">\n";
186 char *tag_End_DataArray = "</DataArray>\n";
187
188 const int VTK_HEX20_INDEX[] =
189 { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 16, 17, 18, 19, 12, 13, 14, 15 };
190 const int VTK_REC9_INDEX[] =
191 { 0, 4, 8, 7, 4, 1, 5, 8, 7, 8, 6, 3, 8, 5, 2, 6 };
192 const int VTK_HEX27_INDEX[] =
193 { 0, 8, 20, 11, 12, 21, 26, 24,
194 8, 1, 9, 20, 21, 13, 22, 26,
195 11, 20, 10, 3, 24, 26, 23, 15,
196 20, 9, 2, 10, 26, 22, 14, 23,
197 12, 21, 26, 24, 4, 16, 25, 19,
198 21, 13, 22, 26, 16, 5, 17, 25,
199 24, 26, 23, 15, 19, 25, 18, 7,
200 26, 22, 14, 23, 25, 17, 6, 18 };
201
202 /* if there is no mesh we just return */
203 if (mesh_p==NULL) return;
204
205 nDim = mesh_p->Nodes->numDim;
206
207 if (nDim != 2 && nDim != 3) {
208 Finley_setError(TYPE_ERROR, "saveVTK: spatial dimension 2 or 3 is supported only.");
209 return;
210 }
211 my_mpi_rank = mesh_p->Nodes->MPIInfo->rank;
212 mpi_size = mesh_p->Nodes->MPIInfo->size;
213
214 /************************************************************************
215 * open the file and check handle *
216 */
217 if (mpi_size > 1) {
218 #ifdef PASO_MPI
219 const int amode = MPI_MODE_CREATE|MPI_MODE_WRONLY|MPI_MODE_UNIQUE_OPEN;
220 int ierr;
221 if (my_mpi_rank == 0 && Paso_fileExists(filename_p)) {
222 remove(filename_p);
223 }
224 ierr = MPI_File_open(mesh_p->Nodes->MPIInfo->comm, (char*)filename_p,
225 amode, mpi_info, &mpi_fileHandle_p);
226 if (ierr != MPI_SUCCESS) {
227 sprintf(errorMsg, "saveVTK: File %s could not be opened for writing in parallel.", filename_p);
228 Finley_setError(IO_ERROR, errorMsg);
229 } else {
230 ierr=MPI_File_set_view(mpi_fileHandle_p,MPI_DISPLACEMENT_CURRENT,
231 MPI_CHAR, MPI_CHAR, "native", mpi_info);
232 }
233 #endif /* PASO_MPI */
234 } else {
235 fileHandle_p = fopen(filename_p, "w");
236 if (fileHandle_p==NULL) {
237 sprintf(errorMsg, "saveVTK: File %s could not be opened for writing.", filename_p);
238 Finley_setError(IO_ERROR, errorMsg);
239 }
240 }
241 if (!Paso_MPIInfo_noError(mesh_p->Nodes->MPIInfo)) return;
242
243 /* General note: From this point if an error occurs Finley_setError is
244 * called and subsequent steps are skipped until the end of this function
245 * where allocated memory is freed and the file is closed. */
246
247 /************************************************************************/
248 /* find the mesh type to be written */
249
250 isCellCentered = TMPMEMALLOC(num_data, bool_t);
251 maxNameLen = 0;
252 if (!Finley_checkPtr(isCellCentered)) {
253 for (dataIdx=0; dataIdx<num_data; ++dataIdx) {
254 if (! isEmpty(data_pp[dataIdx])) {
255 switch(getFunctionSpaceType(data_pp[dataIdx]) ) {
256 case FINLEY_NODES:
257 nodeType = (nodeType == FINLEY_REDUCED_NODES) ? FINLEY_REDUCED_NODES : FINLEY_NODES;
258 isCellCentered[dataIdx] = FALSE;
259 if (elementType==FINLEY_UNKNOWN || elementType==FINLEY_ELEMENTS) {
260 elementType = FINLEY_ELEMENTS;
261 } else {
262 Finley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
263 }
264 break;
265 case FINLEY_REDUCED_NODES:
266 nodeType = FINLEY_REDUCED_NODES;
267 isCellCentered[dataIdx] = FALSE;
268 if (elementType==FINLEY_UNKNOWN || elementType==FINLEY_ELEMENTS) {
269 elementType = FINLEY_ELEMENTS;
270 } else {
271 Finley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
272 }
273 break;
274 case FINLEY_REDUCED_ELEMENTS:
275 hasReducedElements = TRUE;
276 case FINLEY_ELEMENTS:
277 isCellCentered[dataIdx] = TRUE;
278 if (elementType==FINLEY_UNKNOWN || elementType==FINLEY_ELEMENTS) {
279 elementType = FINLEY_ELEMENTS;
280 } else {
281 Finley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
282 }
283 break;
284 case FINLEY_REDUCED_FACE_ELEMENTS:
285 hasReducedElements = TRUE;
286 case FINLEY_FACE_ELEMENTS:
287 isCellCentered[dataIdx] = TRUE;
288 if (elementType==FINLEY_UNKNOWN || elementType==FINLEY_FACE_ELEMENTS) {
289 elementType = FINLEY_FACE_ELEMENTS;
290 } else {
291 Finley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
292 }
293 break;
294 case FINLEY_POINTS:
295 isCellCentered[dataIdx]=TRUE;
296 if (elementType==FINLEY_UNKNOWN || elementType==FINLEY_POINTS) {
297 elementType = FINLEY_POINTS;
298 } else {
299 Finley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
300 }
301 break;
302 case FINLEY_REDUCED_CONTACT_ELEMENTS_1:
303 hasReducedElements = TRUE;
304 case FINLEY_CONTACT_ELEMENTS_1:
305 isCellCentered[dataIdx] = TRUE;
306 if (elementType==FINLEY_UNKNOWN || elementType==FINLEY_CONTACT_ELEMENTS_1) {
307 elementType = FINLEY_CONTACT_ELEMENTS_1;
308 } else {
309 Finley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
310 }
311 break;
312 case FINLEY_REDUCED_CONTACT_ELEMENTS_2:
313 hasReducedElements = TRUE;
314 case FINLEY_CONTACT_ELEMENTS_2:
315 isCellCentered[dataIdx] = TRUE;
316 if (elementType==FINLEY_UNKNOWN || elementType==FINLEY_CONTACT_ELEMENTS_1) {
317 elementType = FINLEY_CONTACT_ELEMENTS_1;
318 } else {
319 Finley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
320 }
321 break;
322 default:
323 sprintf(errorMsg, "saveVTK: unknown function space type %d",getFunctionSpaceType(data_pp[dataIdx]));
324 Finley_setError(TYPE_ERROR, errorMsg);
325 }
326 if (isCellCentered[dataIdx]) {
327 writeCellData = TRUE;
328 } else {
329 writePointData = TRUE;
330 }
331 maxNameLen = MAX(maxNameLen, strlen(names_p[dataIdx]));
332 }
333 }
334 }
335
336 /************************************************************************/
337 /* select number of points and the mesh component */
338
339 if (Finley_noError()) {
340 if (nodeType == FINLEY_REDUCED_NODES) {
341 myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
342 myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
343 globalNumPoints = Finley_NodeFile_getGlobalNumReducedNodes(mesh_p->Nodes);
344 globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
345 } else {
346 myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
347 myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
348 globalNumPoints = Finley_NodeFile_getGlobalNumNodes(mesh_p->Nodes);
349 globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
350 }
351 myNumPoints = myLastNode - myFirstNode;
352 if (elementType==FINLEY_UNKNOWN) elementType=FINLEY_ELEMENTS;
353 switch(elementType) {
354 case FINLEY_ELEMENTS:
355 elements = mesh_p->Elements;
356 break;
357 case FINLEY_FACE_ELEMENTS:
358 elements = mesh_p->FaceElements;
359 break;
360 case FINLEY_POINTS:
361 elements = mesh_p->Points;
362 break;
363 case FINLEY_CONTACT_ELEMENTS_1:
364 elements = mesh_p->ContactElements;
365 break;
366 }
367 if (elements==NULL) {
368 Finley_setError(SYSTEM_ERROR, "saveVTK: undefined element file");
369 } else {
370 /* map finley element type to VTK element type */
371 numCells = elements->numElements;
372 globalNumCells = Finley_ElementFile_getGlobalNumElements(elements);
373 myNumCells = Finley_ElementFile_getMyNumElements(elements);
374 myFirstCell = Finley_ElementFile_getFirstElement(elements);
375 NN = elements->numNodes;
376 if (nodeType==FINLEY_REDUCED_NODES) {
377 typeId = elements->LinearReferenceElement->Type->TypeId;
378 if (hasReducedElements) {
379 quadNodes_p=elements->LinearReferenceElementReducedOrder->QuadNodes;
380 } else {
381 quadNodes_p=elements->LinearReferenceElement->QuadNodes;
382 }
383 } else {
384 typeId = elements->ReferenceElement->Type->TypeId;
385 if (hasReducedElements)
386 quadNodes_p=elements->ReferenceElementReducedOrder->QuadNodes;
387 else
388 quadNodes_p=elements->ReferenceElement->QuadNodes;
389 }
390 switch (typeId) {
391 case Point1:
392 case Line2Face:
393 case Line3Face:
394 case Point1_Contact:
395 case Line2Face_Contact:
396 case Line3Face_Contact:
397 cellType = VTK_VERTEX;
398 numVTKNodesPerElement = 1;
399 break;
400
401 case Line2:
402 case Tri3Face:
403 case Rec4Face:
404 case Line2_Contact:
405 case Tri3_Contact:
406 case Tri3Face_Contact:
407 case Rec4Face_Contact:
408 cellType = VTK_LINE;
409 numVTKNodesPerElement = 2;
410 break;
411
412 case Tri3:
413 case Tet4Face:
414 case Tet4Face_Contact:
415 cellType = VTK_TRIANGLE;
416 numVTKNodesPerElement = 3;
417 break;
418
419 case Rec4:
420 case Hex8Face:
421 case Rec4_Contact:
422 case Hex8Face_Contact:
423 cellType = VTK_QUAD;
424 numVTKNodesPerElement = 4;
425 break;
426
427 case Rec9:
428 numCellFactor = 4;
429 cellType = VTK_QUAD;
430 numVTKNodesPerElement = 4;
431 break;
432
433 case Tet4:
434 cellType = VTK_TETRA;
435 numVTKNodesPerElement = 4;
436 break;
437
438 case Hex8:
439 cellType = VTK_HEXAHEDRON;
440 numVTKNodesPerElement = 8;
441 break;
442
443 case Line3:
444 case Tri6Face:
445 case Rec8Face:
446 case Line3_Contact:
447 case Tri6Face_Contact:
448 case Rec8Face_Contact:
449 cellType = VTK_QUADRATIC_EDGE;
450 numVTKNodesPerElement = 3;
451 break;
452
453 case Tri6:
454 case Tet10Face:
455 case Tri6_Contact:
456 case Tet10Face_Contact:
457 cellType = VTK_QUADRATIC_TRIANGLE;
458 numVTKNodesPerElement = 6;
459 break;
460
461 case Rec8:
462 case Hex20Face:
463 case Rec8_Contact:
464 case Hex20Face_Contact:
465 cellType = VTK_QUADRATIC_QUAD;
466 numVTKNodesPerElement = 8;
467 break;
468
469 case Tet10:
470 cellType = VTK_QUADRATIC_TETRA;
471 numVTKNodesPerElement = 10;
472 break;
473
474 case Hex20:
475 cellType = VTK_QUADRATIC_HEXAHEDRON;
476 numVTKNodesPerElement = 20;
477 break;
478
479 case Hex27:
480 numCellFactor = 8;
481 cellType = VTK_HEXAHEDRON;
482 numVTKNodesPerElement = 8;
483 break;
484
485 default:
486 sprintf(errorMsg, "saveVTK: Element type %s is not supported by VTK.", elements->ReferenceElement->Type->Name);
487 Finley_setError(VALUE_ERROR, errorMsg);
488 }
489 }
490 }
491
492 /* allocate enough memory for text buffer */
493
494 txtBufferSize = strlen(vtkHeader) + 3*LEN_INT_FORMAT + (30+3*maxNameLen)+strlen(metadata)+strlen(metadata_schema);
495 if (mpi_size > 1) {
496 txtBufferSize = MAX(txtBufferSize, myNumPoints * LEN_TMP_BUFFER);
497 txtBufferSize = MAX(txtBufferSize, numCellFactor * myNumCells *
498 (LEN_INT_FORMAT * numVTKNodesPerElement + 1));
499 txtBufferSize = MAX(txtBufferSize,
500 numCellFactor * myNumCells * LEN_TENSOR_FORMAT);
501 txtBufferSize = MAX(txtBufferSize, myNumPoints * LEN_TENSOR_FORMAT);
502 }
503 txtBuffer = TMPMEMALLOC(txtBufferSize+1, char);
504
505 /* sets error if memory allocation failed */
506 Finley_checkPtr(txtBuffer);
507
508 /************************************************************************/
509 /* write number of points and the mesh component */
510
511 if (Finley_noError()) {
512 const index_t *nodeIndex;
513 if (FINLEY_REDUCED_NODES == nodeType) {
514 nodeIndex = elements->ReferenceElement->Type->linearNodes;
515 } else if (Rec9 == typeId) {
516 nodeIndex = VTK_REC9_INDEX;
517 } else if (Hex20 == typeId) {
518 nodeIndex = VTK_HEX20_INDEX;
519 } else if (Hex27 == typeId) {
520 nodeIndex = VTK_HEX27_INDEX;
521 } else if (numVTKNodesPerElement != NN) {
522 nodeIndex = elements->ReferenceElement->Type->geoNodes;
523 } else {
524 nodeIndex = NULL;
525 }
526
527 if (strlen(metadata)>0) {
528 if (strlen(metadata_schema)>0) {
529 sprintf(txtBuffer, vtkHeader," ",metadata_schema,metadata,"\n",globalNumPoints, numCellFactor*globalNumCells, 3);
530 } else {
531 sprintf(txtBuffer, vtkHeader,"","",metadata,"\n",globalNumPoints, numCellFactor*globalNumCells, 3);
532 }
533 } else {
534 if (strlen(metadata_schema)>0) {
535 sprintf(txtBuffer, vtkHeader," ",metadata_schema,"","",globalNumPoints, numCellFactor*globalNumCells, 3);
536 } else {
537 sprintf(txtBuffer, vtkHeader,"","","","",globalNumPoints, numCellFactor*globalNumCells, 3);
538 }
539 }
540
541 if (mpi_size > 1) {
542 /* write the nodes */
543 MPI_RANK0_WRITE_SHARED(txtBuffer);
544 txtBuffer[0] = '\0';
545 txtBufferInUse = 0;
546 if (nDim==2) {
547 for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
548 if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
549 sprintf(tmpBuffer, VECTOR_FORMAT,
550 mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
551 mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
552 0.);
553 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
554 }
555 }
556 } else {
557 for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
558 if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
559 sprintf(tmpBuffer, VECTOR_FORMAT,
560 mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
561 mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
562 mesh_p->Nodes->Coordinates[INDEX2(2, i, nDim)]);
563 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
564 }
565 }
566 } /* nDim */
567 MPI_WRITE_ORDERED(txtBuffer);
568
569 /* write the cells */
570 MPI_RANK0_WRITE_SHARED(tags_End_Points_and_Start_Conn);
571 txtBuffer[0] = '\0';
572 txtBufferInUse = 0;
573 if (nodeIndex == NULL) {
574 for (i = 0; i < numCells; i++) {
575 if (elements->Owner[i] == my_mpi_rank) {
576 for (j = 0; j < numVTKNodesPerElement; j++) {
577 sprintf(tmpBuffer, INT_FORMAT, globalNodeIndex[elements->Nodes[INDEX2(j, i, NN)]]);
578 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
579 }
580 __STRCAT(txtBuffer, NEWLINE, txtBufferInUse);
581 }
582 }
583 } else {
584 for (i = 0; i < numCells; i++) {
585 if (elements->Owner[i] == my_mpi_rank) {
586 for (l = 0; l < numCellFactor; l++) {
587 const int* idx=&nodeIndex[l*numVTKNodesPerElement];
588 for (j = 0; j < numVTKNodesPerElement; j++) {
589 sprintf(tmpBuffer, INT_FORMAT, globalNodeIndex[elements->Nodes[INDEX2(idx[j], i, NN)]]);
590 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
591 }
592 __STRCAT(txtBuffer, NEWLINE, txtBufferInUse);
593 }
594 }
595 }
596 } /* nodeIndex */
597 MPI_WRITE_ORDERED(txtBuffer);
598
599 /* write the offsets */
600 MPI_RANK0_WRITE_SHARED(tags_End_Conn_and_Start_Offset);
601 txtBuffer[0] = '\0';
602 txtBufferInUse = 0;
603 for (i = numVTKNodesPerElement*(myFirstCell*numCellFactor+1);
604 i <= (myFirstCell+myNumCells)*numVTKNodesPerElement*numCellFactor; i += numVTKNodesPerElement)
605 {
606 sprintf(tmpBuffer, INT_NEWLINE_FORMAT, i);
607 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
608 }
609 MPI_WRITE_ORDERED(txtBuffer);
610
611 /* write element type */
612 sprintf(tmpBuffer, INT_NEWLINE_FORMAT, cellType);
613 MPI_RANK0_WRITE_SHARED(tags_End_Offset_and_Start_Type);
614 txtBuffer[0] = '\0';
615 txtBufferInUse = 0;
616 for (i = numVTKNodesPerElement*(myFirstCell*numCellFactor+1);
617 i <= (myFirstCell+myNumCells)*numVTKNodesPerElement*numCellFactor; i += numVTKNodesPerElement)
618 {
619 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
620 }
621 MPI_WRITE_ORDERED(txtBuffer);
622 /* finalize cell information */
623 strcpy(txtBuffer, "</DataArray>\n</Cells>\n");
624 MPI_RANK0_WRITE_SHARED(txtBuffer);
625
626 } else { /***** mpi_size == 1 *****/
627
628 /* write the nodes */
629 fputs(txtBuffer, fileHandle_p);
630 if (nDim==2) {
631 for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
632 if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
633 fprintf(fileHandle_p, VECTOR_FORMAT,
634 mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
635 mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
636 0.);
637 }
638 }
639 } else {
640 for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
641 if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
642 fprintf(fileHandle_p, VECTOR_FORMAT,
643 mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
644 mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
645 mesh_p->Nodes->Coordinates[INDEX2(2, i, nDim)]);
646 }
647 }
648 } /* nDim */
649
650 /* write the cells */
651 fputs(tags_End_Points_and_Start_Conn, fileHandle_p);
652 if (nodeIndex == NULL) {
653 for (i = 0; i < numCells; i++) {
654 for (j = 0; j < numVTKNodesPerElement; j++) {
655 fprintf(fileHandle_p, INT_FORMAT, globalNodeIndex[elements->Nodes[INDEX2(j, i, NN)]]);
656 }
657 fprintf(fileHandle_p, NEWLINE);
658 }
659 } else {
660 for (i = 0; i < numCells; i++) {
661 for (l = 0; l < numCellFactor; l++) {
662 const int* idx=&nodeIndex[l*numVTKNodesPerElement];
663 for (j = 0; j < numVTKNodesPerElement; j++) {
664 fprintf(fileHandle_p, INT_FORMAT, globalNodeIndex[elements->Nodes[INDEX2(idx[j], i, NN)]]);
665 }
666 fprintf(fileHandle_p, NEWLINE);
667 }
668 }
669 } /* nodeIndex */
670
671 /* write the offsets */
672 fputs(tags_End_Conn_and_Start_Offset, fileHandle_p);
673 for (i = numVTKNodesPerElement; i <= numCells*numVTKNodesPerElement*numCellFactor; i += numVTKNodesPerElement) {
674 fprintf(fileHandle_p, INT_NEWLINE_FORMAT, i);
675 }
676
677 /* write element type */
678 sprintf(tmpBuffer, INT_NEWLINE_FORMAT, cellType);
679 fputs(tags_End_Offset_and_Start_Type, fileHandle_p);
680 for (i = 0; i < numCells*numCellFactor; i++)
681 fputs(tmpBuffer, fileHandle_p);
682 /* finalize cell information */
683 fputs("</DataArray>\n</Cells>\n", fileHandle_p);
684 } /* mpi_size */
685
686 } /* Finley_noError */
687
688 /************************************************************************/
689 /* write cell data */
690
691 if (writeCellData && Finley_noError()) {
692 bool_t set_scalar=FALSE, set_vector=FALSE, set_tensor=FALSE;
693 /* mark the active data arrays */
694 strcpy(txtBuffer, "<CellData");
695 for (dataIdx = 0; dataIdx < num_data; dataIdx++) {
696 if (!isEmpty(data_pp[dataIdx]) && isCellCentered[dataIdx]) {
697 /* rank == 0 <--> scalar data */
698 /* rank == 1 <--> vector data */
699 /* rank == 2 <--> tensor data */
700 switch (getDataPointRank(data_pp[dataIdx])) {
701 case 0:
702 if (!set_scalar) {
703 strcat(txtBuffer, " Scalars=\"");
704 strcat(txtBuffer, names_p[dataIdx]);
705 strcat(txtBuffer, "\"");
706 set_scalar = TRUE;
707 }
708 break;
709 case 1:
710 if (!set_vector) {
711 strcat(txtBuffer, " Vectors=\"");
712 strcat(txtBuffer, names_p[dataIdx]);
713 strcat(txtBuffer, "\"");
714 set_vector = TRUE;
715 }
716 break;
717 case 2:
718 if (!set_tensor) {
719 strcat(txtBuffer, " Tensors=\"");
720 strcat(txtBuffer, names_p[dataIdx]);
721 strcat(txtBuffer, "\"");
722 set_tensor = TRUE;
723 }
724 break;
725 default:
726 sprintf(errorMsg, "saveVTK: data %s: VTK supports data with rank <= 2 only.", names_p[dataIdx]);
727 Finley_setError(VALUE_ERROR, errorMsg);
728 }
729 }
730 if (!Finley_noError())
731 break;
732 }
733 }
734 /* only continue if no error occurred */
735 if (writeCellData && Finley_noError()) {
736 strcat(txtBuffer, ">\n");
737 if ( mpi_size > 1) {
738 MPI_RANK0_WRITE_SHARED(txtBuffer);
739 } else {
740 fputs(txtBuffer, fileHandle_p);
741 }
742
743 /* write the arrays */
744 for (dataIdx = 0; dataIdx < num_data; dataIdx++) {
745 if (!isEmpty(data_pp[dataIdx]) && isCellCentered[dataIdx]) {
746 dim_t numPointsPerSample=getNumDataPointsPerSample(data_pp[dataIdx]);
747 dim_t rank = getDataPointRank(data_pp[dataIdx]);
748 dim_t nComp = getDataPointSize(data_pp[dataIdx]);
749 dim_t nCompReqd = 1; /* number of components mpi_required */
750 if (rank == 0) {
751 nCompReqd = 1;
752 shape = 0;
753 } else if (rank == 1) {
754 nCompReqd = 3;
755 shape = getDataPointShape(data_pp[dataIdx], 0);
756 if (shape > 3) {
757 Finley_setError(VALUE_ERROR, "saveVTK: rank 1 objects must have 3 components at most.");
758 }
759 } else {
760 nCompReqd = 9;
761 shape = getDataPointShape(data_pp[dataIdx], 0);
762 if (shape > 3 || shape != getDataPointShape(data_pp[dataIdx], 1)) {
763 Finley_setError(VALUE_ERROR, "saveVTK: rank 2 objects of shape 2x2 or 3x3 supported only.");
764 }
765 }
766 /* bail out if an error occurred */
767 if (!Finley_noError())
768 break;
769
770 sprintf(txtBuffer, tag_Float_DataArray, names_p[dataIdx], nCompReqd);
771 if ( mpi_size > 1) {
772 MPI_RANK0_WRITE_SHARED(txtBuffer);
773 } else {
774 fputs(txtBuffer, fileHandle_p);
775 }
776
777 txtBuffer[0] = '\0';
778 txtBufferInUse = 0;
779 for (i=0; i<numCells; i++) {
780 if (elements->Owner[i] == my_mpi_rank) {
781 void* sampleBuffer=allocSampleBuffer(data_pp[dataIdx]);
782 __const double *values = getSampleDataRO(data_pp[dataIdx], i,sampleBuffer);
783 for (l = 0; l < numCellFactor; l++) {
784 double sampleAvg[NCOMP_MAX];
785 dim_t nCompUsed = MIN(nComp, NCOMP_MAX);
786
787 /* average over number of points in the sample */
788 if (isExpanded(data_pp[dataIdx])) {
789 dim_t hits=0, hits_old;
790 for (k=0; k<nCompUsed; k++) sampleAvg[k]=0;
791 for (j=0; j<numPointsPerSample; j++) {
792 hits_old=hits;
793 if (nodeInQuadrant(quadNodes_p, typeId, j, l)) {
794 hits++;
795 for (k=0; k<nCompUsed; k++) {
796 sampleAvg[k] += values[INDEX2(k,j,nComp)];
797 }
798 }
799 }
800 for (k=0; k<nCompUsed; k++)
801 sampleAvg[k] /= MAX(hits, 1);
802 } else {
803 for (k=0; k<nCompUsed; k++)
804 sampleAvg[k] = values[k];
805 } /* isExpanded */
806
807 /* if the number of required components is more than
808 * the number of actual components, pad with zeros
809 */
810 /* probably only need to get shape of first element */
811 if (nCompReqd == 1) {
812 sprintf(tmpBuffer, SCALAR_FORMAT, sampleAvg[0]);
813 } else if (nCompReqd == 3) {
814 if (shape==1) {
815 sprintf(tmpBuffer, VECTOR_FORMAT,
816 sampleAvg[0], 0.f, 0.f);
817 } else if (shape==2) {
818 sprintf(tmpBuffer, VECTOR_FORMAT,
819 sampleAvg[0], sampleAvg[1], 0.f);
820 } else if (shape==3) {
821 sprintf(tmpBuffer, VECTOR_FORMAT,
822 sampleAvg[0],sampleAvg[1],sampleAvg[2]);
823 }
824 } else if (nCompReqd == 9) {
825 if (shape==1) {
826 sprintf(tmpBuffer, TENSOR_FORMAT,
827 sampleAvg[0], 0.f, 0.f,
828 0.f, 0.f, 0.f,
829 0.f, 0.f, 0.f);
830 } else if (shape==2) {
831 sprintf(tmpBuffer, TENSOR_FORMAT,
832 sampleAvg[0], sampleAvg[1], 0.f,
833 sampleAvg[2], sampleAvg[3], 0.f,
834 0.f, 0.f, 0.f);
835 } else if (shape==3) {
836 sprintf(tmpBuffer, TENSOR_FORMAT,
837 sampleAvg[0],sampleAvg[1],sampleAvg[2],
838 sampleAvg[3],sampleAvg[4],sampleAvg[5],
839 sampleAvg[6],sampleAvg[7],sampleAvg[8]);
840 }
841 }
842 if ( mpi_size > 1) {
843 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
844 } else {
845 fputs(tmpBuffer, fileHandle_p);
846 }
847 } /* for l (numCellFactor) */
848 freeSampleBuffer(sampleBuffer);
849 } /* if I am the owner */
850 } /* for i (numCells) */
851
852 if ( mpi_size > 1) {
853 MPI_WRITE_ORDERED(txtBuffer);
854 MPI_RANK0_WRITE_SHARED(tag_End_DataArray);
855 } else {
856 fputs(tag_End_DataArray, fileHandle_p);
857 }
858 } /* !isEmpty && cellCentered */
859 } /* for dataIdx */
860
861 strcpy(txtBuffer, "</CellData>\n");
862 if ( mpi_size > 1) {
863 MPI_RANK0_WRITE_SHARED(txtBuffer);
864 } else {
865 fputs(txtBuffer, fileHandle_p);
866 }
867 } /* if noError && writeCellData */
868
869 /************************************************************************/
870 /* write point data */
871
872 if (writePointData && Finley_noError()) {
873 /* mark the active data arrays */
874 bool_t set_scalar=FALSE, set_vector=FALSE, set_tensor=FALSE;
875 strcpy(txtBuffer, "<PointData");
876 for (dataIdx = 0; dataIdx<num_data; dataIdx++) {
877 if (!isEmpty(data_pp[dataIdx]) && !isCellCentered[dataIdx]) {
878 switch (getDataPointRank(data_pp[dataIdx])) {
879 case 0:
880 if (!set_scalar) {
881 strcat(txtBuffer, " Scalars=\"");
882 strcat(txtBuffer, names_p[dataIdx]);
883 strcat(txtBuffer, "\"");
884 set_scalar = TRUE;
885 }
886 break;
887 case 1:
888 if (!set_vector) {
889 strcat(txtBuffer, " Vectors=\"");
890 strcat(txtBuffer, names_p[dataIdx]);
891 strcat(txtBuffer, "\"");
892 set_vector = TRUE;
893 }
894 break;
895 case 2:
896 if (!set_tensor) {
897 strcat(txtBuffer, " Tensors=\"");
898 strcat(txtBuffer, names_p[dataIdx]);
899 strcat(txtBuffer, "\"");
900 set_tensor = TRUE;
901 }
902 break;
903 default:
904 sprintf(errorMsg, "saveVTK: data %s: VTK supports data with rank <= 2 only.", names_p[dataIdx]);
905 Finley_setError(VALUE_ERROR, errorMsg);
906 }
907 }
908 if (!Finley_noError())
909 break;
910 }
911 }
912 /* only continue if no error occurred */
913 if (writePointData && Finley_noError()) {
914 strcat(txtBuffer, ">\n");
915 if ( mpi_size > 1) {
916 MPI_RANK0_WRITE_SHARED(txtBuffer);
917 } else {
918 fputs(txtBuffer, fileHandle_p);
919 }
920
921 /* write the arrays */
922 for (dataIdx = 0; dataIdx < num_data; dataIdx++) {
923 if (!isEmpty(data_pp[dataIdx]) && !isCellCentered[dataIdx]) {
924 Finley_NodeMapping* nodeMapping;
925 dim_t rank = getDataPointRank(data_pp[dataIdx]);
926 dim_t nCompReqd = 1; /* number of components mpi_required */
927 if (getFunctionSpaceType(data_pp[dataIdx]) == FINLEY_REDUCED_NODES) {
928 nodeMapping = mesh_p->Nodes->reducedNodesMapping;
929 } else {
930 nodeMapping = mesh_p->Nodes->nodesMapping;
931 }
932 if (rank == 0) {
933 nCompReqd = 1;
934 shape = 0;
935 } else if (rank == 1) {
936 nCompReqd = 3;
937 shape = getDataPointShape(data_pp[dataIdx], 0);
938 if (shape > 3) {
939 Finley_setError(VALUE_ERROR, "saveVTK: rank 1 objects must have 3 components at most.");
940 }
941 } else {
942 nCompReqd = 9;
943 shape=getDataPointShape(data_pp[dataIdx], 0);
944 if (shape > 3 || shape != getDataPointShape(data_pp[dataIdx], 1)) {
945 Finley_setError(VALUE_ERROR, "saveVTK: rank 2 objects of shape 2x2 or 3x3 supported only.");
946 }
947 }
948 /* bail out if an error occurred */
949 if (!Finley_noError())
950 break;
951
952 sprintf(txtBuffer, tag_Float_DataArray, names_p[dataIdx], nCompReqd);
953 if ( mpi_size > 1) {
954 MPI_RANK0_WRITE_SHARED(txtBuffer);
955 } else {
956 fputs(txtBuffer, fileHandle_p);
957 }
958
959 txtBuffer[0] = '\0';
960 txtBufferInUse = 0;
961 for (i=0; i<mesh_p->Nodes->numNodes; i++) {
962 k = globalNodeIndex[i];
963 if ( (myFirstNode <= k) && (k < myLastNode) ) {
964 void* sampleBuffer=allocSampleBuffer(data_pp[dataIdx]);
965 __const double *values = getSampleDataRO(data_pp[dataIdx], nodeMapping->target[i], sampleBuffer);
966 /* if the number of mpi_required components is more than
967 * the number of actual components, pad with zeros.
968 * Probably only need to get shape of first element */
969 if (nCompReqd == 1) {
970 sprintf(tmpBuffer, SCALAR_FORMAT, values[0]);
971 } else if (nCompReqd == 3) {
972 if (shape==1) {
973 sprintf(tmpBuffer, VECTOR_FORMAT,
974 values[0], 0.f, 0.f);
975 } else if (shape==2) {
976 sprintf(tmpBuffer, VECTOR_FORMAT,
977 values[0], values[1], 0.f);
978 } else if (shape==3) {
979 sprintf(tmpBuffer, VECTOR_FORMAT,
980 values[0], values[1], values[2]);
981 }
982 } else if (nCompReqd == 9) {
983 if (shape==1) {
984 sprintf(tmpBuffer, TENSOR_FORMAT,
985 values[0], 0.f, 0.f,
986 0.f, 0.f, 0.f,
987 0.f, 0.f, 0.f);
988 } else if (shape==2) {
989 sprintf(tmpBuffer, TENSOR_FORMAT,
990 values[0], values[1], 0.f,
991 values[2], values[3], 0.f,
992 0.f, 0.f, 0.f);
993 } else if (shape==3) {
994 sprintf(tmpBuffer, TENSOR_FORMAT,
995 values[0], values[1], values[2],
996 values[3], values[4], values[5],
997 values[6], values[7], values[8]);
998 }
999 }
1000 if ( mpi_size > 1) {
1001 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
1002 } else {
1003 fputs(tmpBuffer, fileHandle_p);
1004 }
1005 freeSampleBuffer(sampleBuffer); /* no-one needs values anymore */
1006 } /* if this is my node */
1007 } /* for i (numNodes) */
1008
1009 if ( mpi_size > 1) {
1010 MPI_WRITE_ORDERED(txtBuffer);
1011 MPI_RANK0_WRITE_SHARED(tag_End_DataArray);
1012 } else {
1013 fputs(tag_End_DataArray, fileHandle_p);
1014 }
1015 } /* !isEmpty && !isCellCentered */
1016 } /* for dataIdx */
1017
1018 strcpy(txtBuffer, "</PointData>\n");
1019 if ( mpi_size > 1) {
1020 MPI_RANK0_WRITE_SHARED(txtBuffer);
1021 } else {
1022 fputs(txtBuffer, fileHandle_p);
1023 }
1024 } /* if noError && writePointData */
1025
1026 /* Final write to VTK file */
1027 if (Finley_noError()) {
1028 if (mpi_size > 1) {
1029 MPI_RANK0_WRITE_SHARED(vtkFooter);
1030 } else {
1031 fputs(vtkFooter, fileHandle_p);
1032 }
1033 }
1034
1035 if ( mpi_size > 1) {
1036 #ifdef PASO_MPI
1037 MPI_File_close(&mpi_fileHandle_p);
1038 MPI_Barrier(mesh_p->Nodes->MPIInfo->comm);
1039 #endif
1040 } else {
1041 fclose(fileHandle_p);
1042 }
1043 TMPMEMFREE(isCellCentered);
1044 TMPMEMFREE(txtBuffer);
1045 }
1046

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26