/[escript]/trunk/finley/src/Mesh_saveVTK.c
ViewVC logotype

Contents of /trunk/finley/src/Mesh_saveVTK.c

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26