1 |
/* |
2 |
************************************************************ |
3 |
* Copyright 2006 by ACcESS MNRF * |
4 |
* * |
5 |
* http://www.access.edu.au * |
6 |
* Primary Business: Queensland, Australia * |
7 |
* Licensed under the Open Software License version 3.0 * |
8 |
* http://www.opensource.org/licenses/osl-3.0.php * |
9 |
* * |
10 |
************************************************************ |
11 |
*/ |
12 |
|
13 |
|
14 |
/**************************************************************/ |
15 |
|
16 |
/* writes data and mesh in a vtk file */ |
17 |
|
18 |
/**************************************************************/ |
19 |
|
20 |
/* Author: Paul Cochrane, cochrane@esscc.uq.edu.au */ |
21 |
/* MPI-IO version: Derick Hawcroft, d.hawcroft@uq.edu.au */ |
22 |
|
23 |
/* Version: $Id$ */ |
24 |
|
25 |
/**************************************************************/ |
26 |
|
27 |
|
28 |
#include "Mesh.h" |
29 |
#include "vtkCellType.h" /* copied from vtk source directory !!! */ |
30 |
|
31 |
/* |
32 |
MPI version notes: |
33 |
|
34 |
****************************************************************************** |
35 |
*** **** |
36 |
*** WARNING: Won't work for meshes with peridodic boundary conditions yet **** |
37 |
*** **** |
38 |
****************************************************************************** |
39 |
|
40 |
In this version, the rank==0 process writes *all* opening and closing |
41 |
XML tags. |
42 |
Individual process data is copied to a buffer before being written |
43 |
out. The routines are collectively called and will be called in the natural |
44 |
ordering i.e 0 to maxProcs-1. |
45 |
|
46 |
Notable Notables: |
47 |
the struct localIndexCache: stores local domain indices for faster reference |
48 |
*/ |
49 |
|
50 |
#ifdef PASO_MPI |
51 |
|
52 |
|
53 |
//#define MPIO_HINTS |
54 |
|
55 |
|
56 |
|
57 |
#define MPIO_DEBUG(str) \ |
58 |
{ \ |
59 |
if(myRank == 0) \ |
60 |
printf("==== MPI-IO => %s \n", str); \ |
61 |
} |
62 |
|
63 |
void Finley_Mesh_saveVTK_MPIO(const char * filename_p, Finley_Mesh *mesh_p, const dim_t num_data,char* *names_p, escriptDataC* *data_pp) |
64 |
{ |
65 |
int numPoints, |
66 |
numCells = -1, |
67 |
myRank,comm,gsize, |
68 |
numLocal, |
69 |
nDim, |
70 |
shape; |
71 |
size_t __n; |
72 |
int i,j,k,m,n,count; |
73 |
int numGlobalCells = 0; |
74 |
index_t *nodesGlobal=NULL; // used to get the connectivity right for VTK |
75 |
|
76 |
/* variables associatted with write_celldata/pointdata */ |
77 |
int numPointsPerSample, |
78 |
nComp, |
79 |
nCompReqd; |
80 |
double* values, rtmp; |
81 |
|
82 |
// Local element info (for debugging) |
83 |
size_t numLocalCells, |
84 |
numInternalCells, |
85 |
numBoundaryCells; |
86 |
|
87 |
int rank; |
88 |
|
89 |
int amode = MPI_MODE_CREATE | MPI_MODE_WRONLY | MPI_MODE_SEQUENTIAL; |
90 |
|
91 |
comm = mesh_p->Nodes->MPIInfo->comm; |
92 |
myRank = mesh_p->Nodes->MPIInfo->rank; |
93 |
gsize = mesh_p->Nodes->MPIInfo->size; |
94 |
|
95 |
MPI_File fh; |
96 |
MPI_Status status; |
97 |
MPI_Request req; |
98 |
MPI_Info infoHints; |
99 |
|
100 |
char error_msg[LenErrorMsg_MAX]; |
101 |
|
102 |
int i_data; |
103 |
|
104 |
int nodetype=FINLEY_DEGREES_OF_FREEDOM; |
105 |
int elementtype=FINLEY_UNKNOWN; |
106 |
bool_t isCellCentered[num_data],write_celldata=FALSE,write_pointdata=FALSE; |
107 |
|
108 |
ElementTypeId TypeId; |
109 |
|
110 |
int numVTKNodesPerElement; |
111 |
int cellType; |
112 |
char elemTypeStr[32]; |
113 |
|
114 |
Finley_ElementFile* elements=NULL; |
115 |
|
116 |
|
117 |
// Local node info |
118 |
int numInternalNodes, |
119 |
numLocalNodes, |
120 |
numBoundaryNodes, |
121 |
localDOF; // equals to (DOF of Internal Nodes) + (DOF of Boundary Nodes) of local domain |
122 |
|
123 |
|
124 |
nDim = mesh_p->Nodes->numDim; |
125 |
|
126 |
#ifdef MPIO_HINTS |
127 |
MPI_Info_create(&infoHints); |
128 |
// MPI_Info_set(infoHints, "striping_unit", "424288"); |
129 |
// MPI_Info_set(infoHints, "striping_factor", "16"); |
130 |
// MPI_Info_set(infoHints, "collective_buffering", "true"); |
131 |
// MPI_Info_set(infoHints, "cb_block_size", "131072"); |
132 |
// MPI_Info_set(infoHints, "cb_buffer_size", "1048567"); |
133 |
// MPI_Info_set(infoHints, "cb_nodes", "8"); |
134 |
// MPI_Info_set(infoHints, "access_style", "write_once, sequential"); |
135 |
|
136 |
//XFS only |
137 |
// MPI_Info_set(infoHints, "direct_write", "true"); |
138 |
#else |
139 |
infoHints = MPI_INFO_NULL; |
140 |
#endif |
141 |
|
142 |
// Holds a local node/element index into the global array |
143 |
struct localIndexCache |
144 |
{ |
145 |
index_t *values; |
146 |
int size; |
147 |
}; |
148 |
|
149 |
struct localIndexCache nodeCache, |
150 |
elementCache; |
151 |
|
152 |
// Collective Call |
153 |
MPI_File_open(mesh_p->Nodes->MPIInfo->comm, (char*)filename_p, amode,infoHints, &fh); |
154 |
MPI_File_set_view(fh,MPI_DISPLACEMENT_CURRENT,MPI_CHAR, MPI_CHAR, "native" , infoHints); |
155 |
|
156 |
MPIO_DEBUG(" ***** Enter saveVTK ******") |
157 |
|
158 |
for (i_data=0;i_data<num_data;++i_data) |
159 |
{ |
160 |
if (! isEmpty(data_pp[i_data]) ) |
161 |
{ |
162 |
switch(getFunctionSpaceType(data_pp[i_data])) |
163 |
{ |
164 |
case FINLEY_DEGREES_OF_FREEDOM: |
165 |
nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM; |
166 |
if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) |
167 |
{ |
168 |
elementtype=FINLEY_ELEMENTS; |
169 |
} |
170 |
else |
171 |
{ |
172 |
Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file."); |
173 |
return ; |
174 |
} |
175 |
isCellCentered[i_data]=FALSE; |
176 |
break; |
177 |
case FINLEY_REDUCED_DEGREES_OF_FREEDOM: |
178 |
nodetype = FINLEY_REDUCED_DEGREES_OF_FREEDOM; |
179 |
if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) |
180 |
{ |
181 |
elementtype=FINLEY_ELEMENTS; |
182 |
} |
183 |
else |
184 |
{ |
185 |
Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file."); |
186 |
return; |
187 |
} |
188 |
isCellCentered[i_data]=FALSE; |
189 |
break; |
190 |
case FINLEY_NODES: |
191 |
nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM; |
192 |
if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) |
193 |
{ |
194 |
elementtype=FINLEY_ELEMENTS; |
195 |
} |
196 |
else |
197 |
{ |
198 |
Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file."); |
199 |
return; |
200 |
} |
201 |
isCellCentered[i_data]=FALSE; |
202 |
break; |
203 |
case FINLEY_ELEMENTS: |
204 |
nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM; |
205 |
if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) |
206 |
{ |
207 |
elementtype=FINLEY_ELEMENTS; |
208 |
} |
209 |
else |
210 |
{ |
211 |
Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file."); |
212 |
return; |
213 |
} |
214 |
isCellCentered[i_data]=TRUE; |
215 |
break; |
216 |
case FINLEY_FACE_ELEMENTS: |
217 |
nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM; |
218 |
if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_FACE_ELEMENTS) |
219 |
{ |
220 |
elementtype=FINLEY_FACE_ELEMENTS; |
221 |
} |
222 |
else |
223 |
{ |
224 |
Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file."); |
225 |
return; |
226 |
|
227 |
} |
228 |
isCellCentered[i_data]=TRUE; |
229 |
break; |
230 |
case FINLEY_POINTS: |
231 |
nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM; |
232 |
if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_POINTS) |
233 |
{ |
234 |
elementtype=FINLEY_POINTS; |
235 |
} |
236 |
else |
237 |
{ |
238 |
Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file."); |
239 |
return; |
240 |
|
241 |
} |
242 |
isCellCentered[i_data]=TRUE; |
243 |
break; |
244 |
case FINLEY_CONTACT_ELEMENTS_1: |
245 |
nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM; |
246 |
if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) |
247 |
{ |
248 |
elementtype=FINLEY_CONTACT_ELEMENTS_1; |
249 |
} |
250 |
else |
251 |
{ |
252 |
Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file."); |
253 |
return; |
254 |
|
255 |
} |
256 |
isCellCentered[i_data]=TRUE; |
257 |
break; |
258 |
case FINLEY_CONTACT_ELEMENTS_2: |
259 |
nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM; |
260 |
if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) |
261 |
{ |
262 |
elementtype=FINLEY_CONTACT_ELEMENTS_1; |
263 |
} |
264 |
else |
265 |
{ |
266 |
Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file."); |
267 |
return; |
268 |
|
269 |
} |
270 |
isCellCentered[i_data]=TRUE; |
271 |
break; |
272 |
default: |
273 |
sprintf(error_msg,"saveVTK: Finley does not know anything about function space type %d",getFunctionSpaceType(data_pp[i_data])); |
274 |
Finley_setError(TYPE_ERROR,error_msg); |
275 |
return; |
276 |
|
277 |
} |
278 |
|
279 |
if (isCellCentered[i_data]) |
280 |
{ |
281 |
write_celldata=TRUE; |
282 |
} |
283 |
else |
284 |
{ |
285 |
write_pointdata=TRUE; |
286 |
} |
287 |
} |
288 |
} |
289 |
|
290 |
Finley_NodeDistribution *dist; |
291 |
if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) |
292 |
{ |
293 |
dist = mesh_p->Nodes->reducedDegreeOfFreedomDistribution; |
294 |
} |
295 |
else |
296 |
{ |
297 |
dist = mesh_p->Nodes->degreeOfFreedomDistribution; |
298 |
} |
299 |
|
300 |
numInternalNodes = dist->numInternal; |
301 |
numBoundaryNodes = dist->numBoundary; |
302 |
|
303 |
localDOF = dist->numLocal; |
304 |
|
305 |
numPoints = dist->numGlobal; |
306 |
|
307 |
if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS; |
308 |
switch(elementtype) |
309 |
{ |
310 |
case FINLEY_ELEMENTS: |
311 |
elements=mesh_p->Elements; |
312 |
break; |
313 |
case FINLEY_FACE_ELEMENTS: |
314 |
elements=mesh_p->FaceElements; |
315 |
break; |
316 |
case FINLEY_POINTS: |
317 |
elements=mesh_p->Points; |
318 |
break; |
319 |
case FINLEY_CONTACT_ELEMENTS_1: |
320 |
elements=mesh_p->ContactElements; |
321 |
break; |
322 |
} |
323 |
if (elements==NULL) |
324 |
{ |
325 |
Finley_setError(SYSTEM_ERROR,"saveVTK: undefined element file"); |
326 |
return ; |
327 |
} |
328 |
|
329 |
numCells = elements->numElements; |
330 |
numGlobalCells = elements->elementDistribution->vtxdist[gsize]; |
331 |
numLocalCells = elements->elementDistribution->numLocal; |
332 |
numInternalCells = elements->elementDistribution->numInternal; |
333 |
numBoundaryCells = elements->elementDistribution->numBoundary; |
334 |
|
335 |
if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) |
336 |
{ |
337 |
TypeId = elements->LinearReferenceElement->Type->TypeId; |
338 |
} |
339 |
else |
340 |
{ |
341 |
TypeId = elements->ReferenceElement->Type->TypeId; |
342 |
} |
343 |
|
344 |
switch(TypeId) |
345 |
{ |
346 |
case Point1: |
347 |
case Line2Face: |
348 |
case Line3Face: |
349 |
case Point1_Contact: |
350 |
case Line2Face_Contact: |
351 |
case Line3Face_Contact: |
352 |
cellType = VTK_VERTEX; |
353 |
numVTKNodesPerElement = 1; |
354 |
strcpy(elemTypeStr, "VTK_VERTEX"); |
355 |
break; |
356 |
|
357 |
case Line2: |
358 |
case Tri3Face: |
359 |
case Rec4Face: |
360 |
case Line2_Contact: |
361 |
case Tri3_Contact: |
362 |
case Tri3Face_Contact: |
363 |
case Rec4Face_Contact: |
364 |
cellType = VTK_LINE; |
365 |
numVTKNodesPerElement = 2; |
366 |
strcpy(elemTypeStr, "VTK_LINE"); |
367 |
break; |
368 |
|
369 |
case Tri3: |
370 |
case Tet4Face: |
371 |
case Tet4Face_Contact: |
372 |
cellType = VTK_TRIANGLE; |
373 |
numVTKNodesPerElement = 3; |
374 |
strcpy(elemTypeStr, "VTK_TRIANGLE"); |
375 |
break; |
376 |
|
377 |
case Rec4: |
378 |
case Hex8Face: |
379 |
case Rec4_Contact: |
380 |
case Hex8Face_Contact: |
381 |
cellType = VTK_QUAD; |
382 |
numVTKNodesPerElement = 4; |
383 |
strcpy(elemTypeStr, "VTK_QUAD"); |
384 |
break; |
385 |
|
386 |
case Tet4: |
387 |
cellType = VTK_TETRA; |
388 |
numVTKNodesPerElement = 4; |
389 |
strcpy(elemTypeStr, "VTK_TETRA"); |
390 |
break; |
391 |
|
392 |
case Hex8: |
393 |
cellType = VTK_HEXAHEDRON; |
394 |
numVTKNodesPerElement = 8; |
395 |
strcpy(elemTypeStr, "VTK_HEXAHEDRON"); |
396 |
break; |
397 |
|
398 |
case Line3: |
399 |
case Tri6Face: |
400 |
case Rec8Face: |
401 |
case Line3_Contact: |
402 |
case Tri6Face_Contact: |
403 |
case Rec8Face_Contact: |
404 |
cellType = VTK_QUADRATIC_EDGE; |
405 |
numVTKNodesPerElement = 3; |
406 |
strcpy(elemTypeStr, "VTK_QUADRATIC_EDGE"); |
407 |
break; |
408 |
|
409 |
case Tri6: |
410 |
case Tet10Face: |
411 |
case Tri6_Contact: |
412 |
case Tet10Face_Contact: |
413 |
cellType = VTK_QUADRATIC_TRIANGLE; |
414 |
numVTKNodesPerElement = 6; |
415 |
strcpy(elemTypeStr, "VTK_QUADRATIC_TRIANGLE"); |
416 |
break; |
417 |
|
418 |
case Rec8: |
419 |
case Hex20Face: |
420 |
case Rec8_Contact: |
421 |
case Hex20Face_Contact: |
422 |
cellType = VTK_QUADRATIC_QUAD; |
423 |
numVTKNodesPerElement = 8; |
424 |
strcpy(elemTypeStr, "VTK_QUADRATIC_QUAD"); |
425 |
break; |
426 |
|
427 |
case Tet10: |
428 |
cellType = VTK_QUADRATIC_TETRA; |
429 |
numVTKNodesPerElement = 10; |
430 |
strcpy(elemTypeStr, "VTK_QUADRATIC_TETRA"); |
431 |
break; |
432 |
|
433 |
case Hex20: |
434 |
cellType = VTK_QUADRATIC_HEXAHEDRON; |
435 |
numVTKNodesPerElement = 20; |
436 |
strcpy(elemTypeStr, "VTK_QUADRATIC_HEXAHEDRON"); |
437 |
break; |
438 |
|
439 |
default: |
440 |
sprintf(error_msg, "saveVTK: Element type %s is not supported by VTK",elements->ReferenceElement->Type->Name); |
441 |
Finley_setError(VALUE_ERROR,error_msg); |
442 |
return; |
443 |
} |
444 |
|
445 |
/* Write XML Header */ |
446 |
if(myRank == 0) |
447 |
{ |
448 |
char header[400]; |
449 |
|
450 |
sprintf(header,"<?xml version=\"1.0\"?>\n" \ |
451 |
"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\">\n" \ |
452 |
"<UnstructuredGrid>\n" \ |
453 |
"<Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n" \ |
454 |
"<Points>\n" \ |
455 |
"<DataArray NumberOfComponents=\"%d\" type=\"Float32\" format=\"ascii\">\n" |
456 |
,numPoints,numGlobalCells,MAX(3,nDim)); |
457 |
|
458 |
|
459 |
MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req); |
460 |
MPI_Wait(&req,&status); |
461 |
} |
462 |
|
463 |
MPIO_DEBUG(" Writing Coordinate Points... ") |
464 |
|
465 |
numLocalNodes=localDOF; |
466 |
|
467 |
// we will be writing values which vary from 13-15 chars hence the strlen() |
468 |
char* largebuf = MEMALLOC( numLocalNodes*15*nDim + numLocalNodes*2 + 1 ,char); |
469 |
largebuf[0] = '\0'; |
470 |
char tmpbuf[15]; |
471 |
int tsz=0; |
472 |
int numNodesOutput=0; |
473 |
index_t pos=0; |
474 |
|
475 |
index_t *vtxdist = NULL, *DOFNodes=NULL,*forwardBuffer=NULL,*backwardBuffer=NULL; |
476 |
|
477 |
DOFNodes = MEMALLOC(mesh_p->Nodes->numNodes,index_t); |
478 |
nodeCache.values = MEMALLOC( numLocalNodes, index_t); |
479 |
index_t bc_pos = 0; |
480 |
|
481 |
// Custom string concat: avoids expensive strlen(3) call by strcat(3) |
482 |
// Note the implicit assumption on the variable "tsz" |
483 |
int __len,__j; |
484 |
char *zero = "0.000000e+00 "; |
485 |
char *newline = "\n"; |
486 |
|
487 |
#define __STRCAT(dest,chunk,tsz) \ |
488 |
{ \ |
489 |
__len = strlen(chunk); \ |
490 |
__j = -1; \ |
491 |
while(__j++ < __len) \ |
492 |
*(dest+tsz+__j)=*(chunk+__j); \ |
493 |
tsz+=__len; \ |
494 |
} |
495 |
|
496 |
// Loop over all nodes |
497 |
for (i = 0; i < mesh_p->Nodes->numNodes; i++) |
498 |
{ |
499 |
// This is the bit that will break for periodic BCs because it assumes that there is a one to one |
500 |
// correspondance between nodes and Degrees of freedom |
501 |
//TODO: handle periodic BC's |
502 |
DOFNodes[mesh_p->Nodes->degreeOfFreedom[i]] = i; |
503 |
|
504 |
// Is this node local to the domain ? |
505 |
if( mesh_p->Nodes->degreeOfFreedom[i] < localDOF ) |
506 |
{ |
507 |
for (j = 0; j < nDim; j++) |
508 |
{ |
509 |
sprintf(tmpbuf,"%e ", mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)] ); |
510 |
__STRCAT(largebuf,tmpbuf,tsz) |
511 |
} |
512 |
for (k=0; k<3-nDim; k++) |
513 |
{ |
514 |
__STRCAT(largebuf,zero,tsz) |
515 |
} |
516 |
__STRCAT(largebuf,newline,tsz) |
517 |
nodeCache.values[numNodesOutput++]=i; |
518 |
} |
519 |
} |
520 |
|
521 |
nodeCache.size=numNodesOutput; |
522 |
|
523 |
largebuf[tsz] = '\0'; |
524 |
MPI_File_write_ordered(fh, largebuf,tsz, MPI_CHAR, &status); |
525 |
MEMFREE(largebuf); |
526 |
|
527 |
nodesGlobal = MEMALLOC(mesh_p->Nodes->numNodes ,index_t); |
528 |
|
529 |
// form distribution info on who output which nodes |
530 |
vtxdist = MEMALLOC( gsize+1, index_t ); |
531 |
vtxdist[0]=0; |
532 |
MPI_Allgather(&numNodesOutput,1,MPI_INT,vtxdist+1,1,MPI_INT,comm); |
533 |
for( i=0; i<gsize; i++ ) |
534 |
vtxdist[i+1]+=vtxdist[i]; |
535 |
|
536 |
// will not work for periodic boundary conditions |
537 |
// calculate the local nodes file positions |
538 |
pos = 0; |
539 |
for( i=0; i<mesh_p->Nodes->numNodes; i++ ) |
540 |
{ |
541 |
if( mesh_p->Nodes->degreeOfFreedom[i]< localDOF ) |
542 |
{ |
543 |
nodesGlobal[i] = vtxdist[myRank] + pos++; |
544 |
} |
545 |
else |
546 |
nodesGlobal[i] = -1; |
547 |
} |
548 |
|
549 |
// communicate the local Nodes file position to the interested parties |
550 |
// send local info |
551 |
forwardBuffer = MEMALLOC( mesh_p->Nodes->numNodes, index_t ); |
552 |
for( n=0; n < dist->numNeighbours; n++ ) |
553 |
{ |
554 |
if( dist->edges[n]->numForward) |
555 |
{ |
556 |
for( i=0; i < dist->edges[n]->numForward; i++ ) |
557 |
forwardBuffer[i] = nodesGlobal[DOFNodes[dist->edges[n]->indexForward[i] ]]; |
558 |
Paso_CommBuffer_pack( mesh_p->Nodes->CommBuffer, dist->neighbours[n], NULL, forwardBuffer, sizeof(index_t), 0 ); |
559 |
Paso_CommBuffer_send( mesh_p->Nodes->CommBuffer, dist->neighbours[n], sizeof(index_t) ); |
560 |
} |
561 |
} |
562 |
// receive external info |
563 |
backwardBuffer = MEMALLOC( mesh_p->Nodes->numNodes, index_t ); |
564 |
for( n=0; n < dist->numNeighbours; n++ ) |
565 |
{ |
566 |
if( dist->edges[n]->numBackward ) |
567 |
{ |
568 |
Paso_CommBuffer_recv(mesh_p->Nodes->CommBuffer, dist->neighbours[n], sizeof(index_t)); |
569 |
Paso_CommBuffer_unpack(mesh_p->Nodes->CommBuffer, dist->neighbours[n], NULL, backwardBuffer, sizeof(index_t), 0 ); |
570 |
/* TODO: voodoo to handle perdiodic BC's */ |
571 |
for( i=0; i<dist->edges[n]->numBackward; i++ ) |
572 |
nodesGlobal[DOFNodes[dist->edges[n]->indexBackward[i] ]] = backwardBuffer[i]; |
573 |
} |
574 |
} |
575 |
|
576 |
|
577 |
|
578 |
MEMFREE(vtxdist); |
579 |
MEMFREE(DOFNodes); |
580 |
MEMFREE(backwardBuffer); |
581 |
MEMFREE(forwardBuffer); |
582 |
|
583 |
if( myRank == 0) |
584 |
{ |
585 |
char* tags = "</DataArray>\n</Points>\n<Cells>\n<DataArray Name=\"connectivity\" type=\"Int32\" " \ |
586 |
"format=\"ascii\">\n" ; |
587 |
MPI_File_iwrite_shared(fh,tags,strlen(tags),MPI_CHAR,&req); |
588 |
MPI_Wait(&req,&status); |
589 |
} |
590 |
MPIO_DEBUG(" Done Writing Coordinate Points ") |
591 |
|
592 |
/* BEGIN CONNECTIVITY */ |
593 |
|
594 |
int NN = elements->ReferenceElement->Type->numNodes; /* num Nodes holding ref-element */ |
595 |
|
596 |
// Collective |
597 |
MPIO_DEBUG(" Writing Connectivity... ") |
598 |
|
599 |
size_t sz = numLocalCells*6*numVTKNodesPerElement + numLocalCells; |
600 |
largebuf = MEMALLOC(sz,char); |
601 |
largebuf[0] = '\0'; |
602 |
tsz=0; |
603 |
pos = 0; |
604 |
// numCells? |
605 |
elementCache.values = MEMALLOC(numLocalCells,index_t); |
606 |
if (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) |
607 |
{ |
608 |
for (i = 0; i < numCells; i++) |
609 |
{ |
610 |
|
611 |
if (elements->Id[i] >= elements->elementDistribution->vtxdist[i] && elements->Id[i] <= elements->elementDistribution->vtxdist[i+1] - 1 ) |
612 |
{ |
613 |
for (j = 0; j < numVTKNodesPerElement; j++) |
614 |
{ |
615 |
sprintf(tmpbuf,"%d ",nodesGlobal[mesh_p->Nodes->toReduced[elements->Nodes[INDEX2(elements->ReferenceElement->Type->linearNodes[j], i, NN)]]]); |
616 |
__STRCAT(largebuf,tmpbuf,tsz) |
617 |
} |
618 |
__STRCAT(largebuf,newline,tsz) |
619 |
elementCache.values[pos++]=i; |
620 |
} |
621 |
} |
622 |
} |
623 |
else if (VTK_QUADRATIC_HEXAHEDRON==cellType) |
624 |
{ |
625 |
char tmpbuf2[20*20*2]; |
626 |
|
627 |
for (i = 0; i < numCells; i++) |
628 |
{ |
629 |
if( elements->Id[i] >= elements->elementDistribution->vtxdist[myRank] && elements->Id[i] <= elements->elementDistribution->vtxdist[myRank+1]-1) |
630 |
{ |
631 |
sprintf(tmpbuf2,"%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d\n", |
632 |
nodesGlobal[elements->Nodes[INDEX2(0, i, NN)]], |
633 |
nodesGlobal[elements->Nodes[INDEX2(1, i, NN)]], |
634 |
nodesGlobal[elements->Nodes[INDEX2(2, i, NN)]], |
635 |
nodesGlobal[elements->Nodes[INDEX2(3, i, NN)]], |
636 |
nodesGlobal[elements->Nodes[INDEX2(4, i, NN)]], |
637 |
nodesGlobal[elements->Nodes[INDEX2(5, i, NN)]], |
638 |
nodesGlobal[elements->Nodes[INDEX2(6, i, NN)]], |
639 |
nodesGlobal[elements->Nodes[INDEX2(7, i, NN)]], |
640 |
nodesGlobal[elements->Nodes[INDEX2(8, i, NN)]], |
641 |
nodesGlobal[elements->Nodes[INDEX2(9, i, NN)]], |
642 |
nodesGlobal[elements->Nodes[INDEX2(10, i, NN)]], |
643 |
nodesGlobal[elements->Nodes[INDEX2(11, i, NN)]], |
644 |
nodesGlobal[elements->Nodes[INDEX2(16, i, NN)]], |
645 |
nodesGlobal[elements->Nodes[INDEX2(17, i, NN)]], |
646 |
nodesGlobal[elements->Nodes[INDEX2(18, i, NN)]], |
647 |
nodesGlobal[elements->Nodes[INDEX2(19, i, NN)]], |
648 |
nodesGlobal[elements->Nodes[INDEX2(12, i, NN)]], |
649 |
nodesGlobal[elements->Nodes[INDEX2(13, i, NN)]], |
650 |
nodesGlobal[elements->Nodes[INDEX2(14, i, NN)]], |
651 |
nodesGlobal[elements->Nodes[INDEX2(15, i, NN)]]); |
652 |
__STRCAT(largebuf,tmpbuf2,tsz) |
653 |
elementCache.values[pos++]=i; |
654 |
} |
655 |
} |
656 |
} |
657 |
else if (numVTKNodesPerElement!=NN) |
658 |
{ |
659 |
|
660 |
for (i = 0; i < numCells; i++) |
661 |
{ |
662 |
if (elements->Id[i] >= elements->elementDistribution->vtxdist[i] && elements->Id[i] <= elements->elementDistribution->vtxdist[i+1] - 1 ) |
663 |
{ |
664 |
for (j = 0; j < numVTKNodesPerElement; j++) |
665 |
{ |
666 |
sprintf(tmpbuf,"%d ", nodesGlobal[elements->Nodes[INDEX2(elements->ReferenceElement->Type->geoNodes[j], i, NN)]]); |
667 |
__STRCAT(largebuf,tmpbuf,tsz) |
668 |
} |
669 |
__STRCAT(largebuf,newline,tsz) |
670 |
elementCache.values[pos++]=i; |
671 |
} |
672 |
} |
673 |
} |
674 |
else |
675 |
{ |
676 |
for(i = 0;i < numCells ; i++) |
677 |
{ |
678 |
if( elements->Id[i] >= elements->elementDistribution->vtxdist[myRank] && elements->Id[i] <= elements->elementDistribution->vtxdist[myRank+1]-1) |
679 |
{ |
680 |
for (j = 0; j < numVTKNodesPerElement; j++) |
681 |
{ |
682 |
sprintf(tmpbuf,"%d ", nodesGlobal[ elements->Nodes[INDEX2(j, i, NN) ] ] ); |
683 |
__STRCAT(largebuf,tmpbuf,tsz) |
684 |
} |
685 |
__STRCAT(largebuf,newline,tsz) |
686 |
elementCache.values[pos++]=i; |
687 |
} |
688 |
} |
689 |
} |
690 |
|
691 |
elementCache.size = pos; |
692 |
|
693 |
largebuf[tsz] = '\0'; |
694 |
MPI_File_write_ordered(fh,largebuf,tsz, MPI_CHAR, &status); |
695 |
MEMFREE(largebuf); |
696 |
MPIO_DEBUG(" Done Writing Connectivity ") |
697 |
MPIO_DEBUG(" Writing Offsets & Types... ") |
698 |
|
699 |
// Non-Collective |
700 |
if( myRank == 0) |
701 |
{ |
702 |
// write out the DataArray element for the offsets |
703 |
char* tag1 = "</DataArray>\n<DataArray Name=\"offsets\" type=\"Int32\" format=\"ascii\">\n"; |
704 |
char* tag2 = "</DataArray>\n"; |
705 |
char *tag3 = "<DataArray Name=\"types\" type=\"UInt8\" format=\"ascii\">\n"; |
706 |
char *tag4 = "</DataArray>\n</Cells>\n"; |
707 |
|
708 |
int n = numVTKNodesPerElement; |
709 |
|
710 |
// allocate an upper bound on number of bytes needed |
711 |
int sz=0; |
712 |
int lg = log10(numGlobalCells * n) + 1; |
713 |
sz += numGlobalCells*lg; |
714 |
sz += numGlobalCells; |
715 |
tsz = 0; |
716 |
|
717 |
char* largebuf = MEMALLOC(sz + strlen(tag1) + strlen(tag2) + strlen(tag3) + strlen(tag4),char); |
718 |
largebuf[0] ='\0'; |
719 |
char tmp[10]; |
720 |
__STRCAT(largebuf,tag1,tsz) |
721 |
for (i=numVTKNodesPerElement; i<=numGlobalCells*numVTKNodesPerElement; i+=numVTKNodesPerElement) |
722 |
{ |
723 |
sprintf(tmp,"%d\n", i); |
724 |
__STRCAT(largebuf,tmp,tsz) |
725 |
} |
726 |
__STRCAT(largebuf,tag2,tsz) |
727 |
largebuf[tsz] = '\0'; |
728 |
MPI_File_iwrite_shared(fh,largebuf, tsz,MPI_CHAR,&req); |
729 |
MPI_Wait(&req,&status); |
730 |
|
731 |
// re-using buffer! |
732 |
largebuf[0] = '\0'; |
733 |
tsz = 0; |
734 |
__STRCAT(largebuf,tag3,tsz) |
735 |
for (i=0; i<numGlobalCells; i++) |
736 |
{ |
737 |
sprintf(tmp, "%d\n", cellType); |
738 |
__STRCAT(largebuf,tmp,tsz) |
739 |
} |
740 |
__STRCAT(largebuf,tag4,tsz) |
741 |
largebuf[tsz] = '\0'; |
742 |
MPI_File_iwrite_shared(fh,largebuf,tsz,MPI_CHAR,&req); |
743 |
MPI_Wait(&req,&status); |
744 |
MEMFREE(largebuf); |
745 |
} |
746 |
|
747 |
MPIO_DEBUG(" Done Writing Offsets & Types ") |
748 |
|
749 |
// Write Cell data header Tags |
750 |
if(myRank == 0) |
751 |
{ |
752 |
MPIO_DEBUG(" Writing Cell Data ...") |
753 |
if( write_celldata) |
754 |
{ |
755 |
char tmpBuf[80]; |
756 |
char header[600]; |
757 |
// mark the active data arrays |
758 |
bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE; |
759 |
sprintf(tmpBuf, "<CellData"); |
760 |
strcat(header,tmpBuf); |
761 |
for (i_data =0 ;i_data<num_data;++i_data) |
762 |
{ |
763 |
if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data]) |
764 |
{ |
765 |
// if the rank == 0: --> scalar data |
766 |
// if the rank == 1: --> vector data |
767 |
// if the rank == 2: --> tensor data |
768 |
|
769 |
switch(getDataPointRank(data_pp[i_data])) |
770 |
{ |
771 |
case 0: |
772 |
if (! set_scalar) |
773 |
{ |
774 |
sprintf(tmpBuf," Scalars=\"%s\"",names_p[i_data]); |
775 |
strcat(header,tmpBuf); |
776 |
set_scalar=TRUE; |
777 |
} |
778 |
break; |
779 |
case 1: |
780 |
if (! set_vector) |
781 |
{ |
782 |
sprintf(tmpBuf," Vectors=\"%s\"",names_p[i_data]); |
783 |
strcat(header,tmpBuf); |
784 |
set_vector=TRUE; |
785 |
} |
786 |
break; |
787 |
case 2: |
788 |
if (! set_tensor) |
789 |
{ |
790 |
sprintf(tmpBuf," Tensors=\"%s\"",names_p[i_data]); |
791 |
strcat(header,tmpBuf); |
792 |
set_tensor=TRUE; |
793 |
} |
794 |
break; |
795 |
default: |
796 |
sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]); |
797 |
Finley_setError(VALUE_ERROR,error_msg); |
798 |
return; |
799 |
} |
800 |
} |
801 |
} |
802 |
strcat(header, ">\n"); |
803 |
MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req); |
804 |
MPI_Wait(&req,&status); |
805 |
} |
806 |
} |
807 |
|
808 |
// write actual data (collective) |
809 |
if(write_celldata) |
810 |
{ |
811 |
for (i_data =0 ;i_data<num_data;++i_data) |
812 |
{ |
813 |
if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data]) |
814 |
{ |
815 |
numPointsPerSample = elements->ReferenceElement->numQuadNodes; |
816 |
rank = getDataPointRank(data_pp[i_data]); |
817 |
nComp = getDataPointSize(data_pp[i_data]); |
818 |
nCompReqd=1; // the number of components required by vtk |
819 |
shape=0; |
820 |
if (rank == 0) |
821 |
{ |
822 |
nCompReqd = 1; |
823 |
} |
824 |
else if (rank == 1) |
825 |
{ |
826 |
shape=getDataPointShape(data_pp[i_data], 0); |
827 |
if (shape>3) |
828 |
{ |
829 |
Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components"); |
830 |
return; |
831 |
} |
832 |
nCompReqd = 3; |
833 |
} |
834 |
else |
835 |
{ |
836 |
shape=getDataPointShape(data_pp[i_data], 0); |
837 |
if (shape>3 || shape != getDataPointShape(data_pp[i_data], 1)) |
838 |
{ |
839 |
Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape"); |
840 |
return; |
841 |
} |
842 |
nCompReqd = 9; |
843 |
} |
844 |
|
845 |
if( myRank == 0) |
846 |
{ |
847 |
char header[250]; |
848 |
sprintf(header,"<DataArray Name=\"%s\" type=\"Float32\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd); |
849 |
MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req); |
850 |
MPI_Wait(&req,&status); |
851 |
} |
852 |
|
853 |
// Write the actual data |
854 |
char tmpbuf[15]; |
855 |
char* largebuf = MEMALLOC(nCompReqd*numLocalCells*15 + numLocalCells*nCompReqd + nCompReqd + 15,char); |
856 |
largebuf[0] = '\0'; |
857 |
size_t tsz = 0; |
858 |
|
859 |
double sampleAvg[nComp]; |
860 |
|
861 |
for (k=0; k<elementCache.size; k++) |
862 |
{ |
863 |
i = elementCache.values[k]; |
864 |
|
865 |
values = getSampleData(data_pp[i_data], i); |
866 |
// averaging over the number of points in the sample |
867 |
for (n=0; n<nComp; n++) |
868 |
{ |
869 |
rtmp = 0.; |
870 |
for (j=0; j<numPointsPerSample; j++) rtmp += values[INDEX2(n,j,nComp)]; |
871 |
sampleAvg[k] = rtmp/numPointsPerSample; |
872 |
} |
873 |
// if the number of required components is more than the number |
874 |
// of actual components, pad with zeros |
875 |
|
876 |
// probably only need to get shape of first element |
877 |
// write the data different ways for scalar, vector and tensor |
878 |
if (nCompReqd == 1) |
879 |
{ |
880 |
sprintf(tmpbuf, " %e", sampleAvg[0]); |
881 |
__STRCAT(largebuf,tmpbuf,tsz) |
882 |
} |
883 |
else if (nCompReqd == 3) |
884 |
{ |
885 |
// write out the data |
886 |
for (m=0; m<shape; m++) |
887 |
{ |
888 |
sprintf(tmpbuf, " %e", sampleAvg[m]); |
889 |
__STRCAT(largebuf,tmpbuf,tsz) |
890 |
} |
891 |
for (m=0; m<nCompReqd-shape; m++) |
892 |
{ |
893 |
__STRCAT(largebuf,zero,tsz) |
894 |
} |
895 |
} |
896 |
else if (nCompReqd == 9) |
897 |
{ |
898 |
// tensor data, so have a 3x3 matrix to output as a row |
899 |
// of 9 data points |
900 |
count = 0; |
901 |
for (m=0; m<shape; m++) |
902 |
{ |
903 |
for (n=0; n<shape; n++) |
904 |
{ |
905 |
sprintf(tmpbuf, " %e", sampleAvg[count]); |
906 |
__STRCAT(largebuf,tmpbuf,tsz) |
907 |
count++; |
908 |
} |
909 |
for (n=0; n<3-shape; n++) |
910 |
{ |
911 |
__STRCAT(largebuf,zero,tsz) |
912 |
} |
913 |
} |
914 |
for (m=0; m<3-shape; m++) |
915 |
for (n=0; n<3; n++) |
916 |
{ |
917 |
__STRCAT(largebuf,zero,tsz) |
918 |
} |
919 |
} |
920 |
__STRCAT(largebuf,newline,tsz) |
921 |
} |
922 |
largebuf[tsz] = '\0'; |
923 |
MPI_File_write_ordered(fh,largebuf,tsz,MPI_CHAR,&status); |
924 |
MEMFREE(largebuf); |
925 |
if( myRank == 0) |
926 |
{ |
927 |
char *tag = "</DataArray>\n"; |
928 |
MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req); |
929 |
MPI_Wait(&req,&status); |
930 |
} |
931 |
|
932 |
} |
933 |
} |
934 |
// closing celldata tag |
935 |
if(myRank == 0) |
936 |
{ |
937 |
char* tag = "</CellData>\n"; |
938 |
MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req); |
939 |
MPI_Wait(&req,&status); |
940 |
} |
941 |
|
942 |
MPIO_DEBUG(" Done Writing Cell Data ") |
943 |
} |
944 |
|
945 |
|
946 |
// Write Point Data Header Tags |
947 |
if( myRank == 0) |
948 |
{ |
949 |
char header[600]; |
950 |
char tmpBuf[50]; |
951 |
|
952 |
if (write_pointdata) |
953 |
{ |
954 |
MPIO_DEBUG(" Writing Pointdata... ") |
955 |
// mark the active data arrays |
956 |
bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE; |
957 |
sprintf(header, "<PointData"); |
958 |
for (i_data =0 ;i_data<num_data;++i_data) |
959 |
{ |
960 |
if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data]) |
961 |
{ |
962 |
// if the rank == 0: --> scalar data |
963 |
// if the rank == 1: --> vector data |
964 |
// if the rank == 2: --> tensor data |
965 |
|
966 |
switch(getDataPointRank(data_pp[i_data])) |
967 |
{ |
968 |
case 0: |
969 |
if (! set_scalar) |
970 |
{ |
971 |
sprintf(tmpBuf," Scalars=\"%s\"",names_p[i_data]); |
972 |
strcat(header,tmpBuf); |
973 |
set_scalar=TRUE; |
974 |
} |
975 |
break; |
976 |
case 1: |
977 |
if (! set_vector) |
978 |
{ |
979 |
sprintf(tmpBuf," Vectors=\"%s\"",names_p[i_data]); |
980 |
strcat(header,tmpBuf); |
981 |
set_vector=TRUE; |
982 |
} |
983 |
break; |
984 |
case 2: |
985 |
if (! set_tensor) |
986 |
{ |
987 |
sprintf(tmpBuf," Tensors=\"%s\"",names_p[i_data]); |
988 |
strcat(header,tmpBuf); |
989 |
set_tensor=TRUE; |
990 |
} |
991 |
break; |
992 |
default: |
993 |
sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]); |
994 |
Finley_setError(VALUE_ERROR,error_msg); |
995 |
return; |
996 |
} |
997 |
} |
998 |
} |
999 |
strcat(header, ">\n"); |
1000 |
MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req); |
1001 |
MPI_Wait(&req,&status); |
1002 |
} |
1003 |
} |
1004 |
|
1005 |
// write actual data |
1006 |
if(write_pointdata) |
1007 |
{ |
1008 |
for (i_data =0 ;i_data<num_data;++i_data) |
1009 |
{ |
1010 |
if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data]) |
1011 |
{ |
1012 |
numPointsPerSample = elements->ReferenceElement->numQuadNodes; |
1013 |
rank = getDataPointRank(data_pp[i_data]); |
1014 |
nComp = getDataPointSize(data_pp[i_data]); |
1015 |
nCompReqd=1; // the number of components required by vtk |
1016 |
shape=0; |
1017 |
if (rank == 0) |
1018 |
{ |
1019 |
nCompReqd = 1; |
1020 |
} |
1021 |
else if (rank == 1) |
1022 |
{ |
1023 |
shape=getDataPointShape(data_pp[i_data], 0); |
1024 |
if (shape>3) |
1025 |
{ |
1026 |
Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components"); |
1027 |
return; |
1028 |
} |
1029 |
nCompReqd = 3; |
1030 |
} |
1031 |
else |
1032 |
{ |
1033 |
shape=getDataPointShape(data_pp[i_data], 0); |
1034 |
if (shape>3 || shape != getDataPointShape(data_pp[i_data], 1)) |
1035 |
{ |
1036 |
Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape"); |
1037 |
return; |
1038 |
} |
1039 |
nCompReqd = 9; |
1040 |
} |
1041 |
|
1042 |
if( myRank == 0) |
1043 |
{ |
1044 |
char header[250]; |
1045 |
sprintf(header, "<DataArray Name=\"%s\" type=\"Float32\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd); |
1046 |
MPI_File_iwrite_shared(fh,header,strlen(header),MPI_CHAR,&req); |
1047 |
MPI_Wait(&req,&status); |
1048 |
} |
1049 |
// write out the data |
1050 |
// if the number of required components is more than the number |
1051 |
// of actual components, pad with zeros |
1052 |
|
1053 |
char tmpbuf[15]; |
1054 |
char* largebuf = MEMALLOC(nCompReqd*numLocalNodes*15 + numLocal*nCompReqd + nCompReqd + 15,char); |
1055 |
largebuf[0] = '\0'; |
1056 |
bool_t do_write=TRUE; |
1057 |
size_t tsz = 0; |
1058 |
|
1059 |
for(k=0;k < nodeCache.size;k++) |
1060 |
{ |
1061 |
i = nodeCache.values[k]; |
1062 |
|
1063 |
if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) |
1064 |
{ |
1065 |
if (mesh_p->Nodes->toReduced[i]>=0) |
1066 |
{ |
1067 |
switch(getFunctionSpaceType(data_pp[i_data])) |
1068 |
{ |
1069 |
case FINLEY_DEGREES_OF_FREEDOM: |
1070 |
values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]); |
1071 |
break; |
1072 |
case FINLEY_REDUCED_DEGREES_OF_FREEDOM: |
1073 |
values = getSampleData(data_pp[i_data],mesh_p->Nodes->reducedDegreeOfFreedom[i]); |
1074 |
break; |
1075 |
case FINLEY_NODES: |
1076 |
values = getSampleData(data_pp[i_data],i); |
1077 |
break; |
1078 |
} |
1079 |
do_write=TRUE; |
1080 |
} |
1081 |
else |
1082 |
{ |
1083 |
do_write=FALSE; |
1084 |
} |
1085 |
} |
1086 |
else |
1087 |
{ |
1088 |
do_write=TRUE; |
1089 |
switch(getFunctionSpaceType(data_pp[i_data])) |
1090 |
{ |
1091 |
case FINLEY_DEGREES_OF_FREEDOM: |
1092 |
values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]); |
1093 |
break; |
1094 |
case FINLEY_NODES: |
1095 |
values = getSampleData(data_pp[i_data],i); |
1096 |
break; |
1097 |
} |
1098 |
} |
1099 |
// write the data different ways for scalar, vector and tensor |
1100 |
if (do_write) |
1101 |
{ |
1102 |
if (nCompReqd == 1) |
1103 |
{ |
1104 |
sprintf(tmpbuf," %e", values[0]); |
1105 |
__STRCAT(largebuf,tmpbuf,tsz) |
1106 |
} |
1107 |
else if (nCompReqd == 3) |
1108 |
{ |
1109 |
for (m=0; m<shape; m++) |
1110 |
{ |
1111 |
|
1112 |
sprintf(tmpbuf," %e",values[m]); |
1113 |
__STRCAT(largebuf,tmpbuf,tsz) |
1114 |
} |
1115 |
for (m=0; m<nCompReqd-shape; m++) |
1116 |
{ |
1117 |
__STRCAT(largebuf,zero,tsz) |
1118 |
} |
1119 |
} |
1120 |
else if (nCompReqd == 9) |
1121 |
{ |
1122 |
// tensor data, so have a 3x3 matrix to output as a row |
1123 |
// of 9 data points |
1124 |
count = 0; |
1125 |
for (m=0; m<shape; m++) |
1126 |
{ |
1127 |
for (n=0; n<shape; n++) |
1128 |
{ |
1129 |
sprintf(tmpbuf," %e", values[count]); |
1130 |
__STRCAT(largebuf,tmpbuf,tsz) |
1131 |
count++; |
1132 |
} |
1133 |
for (n=0; n<3-shape; n++) |
1134 |
{ |
1135 |
__STRCAT(largebuf,zero,tsz) |
1136 |
} |
1137 |
} |
1138 |
for (m=0; m<3-shape; m++) |
1139 |
{ |
1140 |
for (n=0; n<3; n++) |
1141 |
{ |
1142 |
__STRCAT(largebuf,zero,tsz) |
1143 |
} |
1144 |
} |
1145 |
} |
1146 |
__STRCAT(largebuf,newline,tsz) |
1147 |
} |
1148 |
|
1149 |
} |
1150 |
// Write out local data |
1151 |
|
1152 |
largebuf[tsz] = '\0'; |
1153 |
MPI_File_write_ordered(fh,largebuf,tsz,MPI_CHAR,&status); |
1154 |
MEMFREE(largebuf); |
1155 |
if( myRank == 0) |
1156 |
{ |
1157 |
char *tag = "</DataArray>\n"; |
1158 |
MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req); |
1159 |
MPI_Wait(&req,&status); |
1160 |
} |
1161 |
} |
1162 |
} |
1163 |
// Finish off with closing tag |
1164 |
if(myRank == 0) |
1165 |
{ |
1166 |
char* tag = "</PointData>\n"; |
1167 |
MPI_File_iwrite_shared(fh,tag,strlen(tag),MPI_CHAR,&req); |
1168 |
MPI_Wait(&req,&status); |
1169 |
} |
1170 |
MPIO_DEBUG(" Done Writing Pointdata ") |
1171 |
} |
1172 |
// end write_pointdata |
1173 |
|
1174 |
// tag and bag... |
1175 |
if (myRank == 0) |
1176 |
{ |
1177 |
char *footer = "</Piece>\n</UnstructuredGrid>\n</VTKFile>"; |
1178 |
MPI_File_iwrite_shared(fh,footer,strlen(footer),MPI_CHAR,&req); |
1179 |
MPI_Wait(&req,&status); |
1180 |
} |
1181 |
|
1182 |
MEMFREE(nodesGlobal); |
1183 |
MEMFREE(nodeCache.values); |
1184 |
MEMFREE(elementCache.values); |
1185 |
#ifdef MPIO_HINTS |
1186 |
MPI_Info_free(&infoHints); |
1187 |
#undef MPIO_HINTS |
1188 |
#endif |
1189 |
MPI_File_close(&fh); |
1190 |
MPIO_DEBUG(" ***** Exit saveVTK ***** ") |
1191 |
#undef __STRCAT |
1192 |
} |
1193 |
|
1194 |
#undef MPIO_DEBUG |
1195 |
#else |
1196 |
|
1197 |
|
1198 |
|
1199 |
|
1200 |
void Finley_Mesh_saveVTK(const char * filename_p, Finley_Mesh *mesh_p, const dim_t num_data,char* *names_p, escriptDataC* *data_pp) |
1201 |
{ |
1202 |
char error_msg[LenErrorMsg_MAX]; |
1203 |
/* if there is no mesh we just return */ |
1204 |
if (mesh_p==NULL) return; |
1205 |
|
1206 |
int i, j, k, numVTKNodesPerElement,i_data,m, count, n, rank,shape, numPoints, cellType, numCells, |
1207 |
nDim, numPointsPerSample, nComp, nCompReqd; |
1208 |
|
1209 |
index_t j2; |
1210 |
double* values, rtmp; |
1211 |
char elemTypeStr[32]; |
1212 |
|
1213 |
/* open the file and check handle */ |
1214 |
|
1215 |
FILE * fileHandle_p = fopen(filename_p, "w"); |
1216 |
if (fileHandle_p==NULL) |
1217 |
{ |
1218 |
sprintf(error_msg, "saveVTK: File %s could not be opened for writing.", filename_p); |
1219 |
Finley_setError(IO_ERROR,error_msg); |
1220 |
return; |
1221 |
} |
1222 |
/* find the mesh type to be written */ |
1223 |
int nodetype=FINLEY_DEGREES_OF_FREEDOM; |
1224 |
int elementtype=FINLEY_UNKNOWN; |
1225 |
bool_t isCellCentered[num_data],write_celldata=FALSE,write_pointdata=FALSE; |
1226 |
for (i_data=0;i_data<num_data;++i_data) |
1227 |
{ |
1228 |
if (! isEmpty(data_pp[i_data])) |
1229 |
{ |
1230 |
switch(getFunctionSpaceType(data_pp[i_data])) |
1231 |
{ |
1232 |
case FINLEY_DEGREES_OF_FREEDOM: |
1233 |
nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM; |
1234 |
if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) |
1235 |
{ |
1236 |
elementtype=FINLEY_ELEMENTS; |
1237 |
} |
1238 |
else |
1239 |
{ |
1240 |
Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file."); |
1241 |
fclose(fileHandle_p); |
1242 |
return; |
1243 |
} |
1244 |
isCellCentered[i_data]=FALSE; |
1245 |
break; |
1246 |
case FINLEY_REDUCED_DEGREES_OF_FREEDOM: |
1247 |
nodetype = FINLEY_REDUCED_DEGREES_OF_FREEDOM; |
1248 |
if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) |
1249 |
{ |
1250 |
elementtype=FINLEY_ELEMENTS; |
1251 |
} |
1252 |
else |
1253 |
{ |
1254 |
Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file."); |
1255 |
fclose(fileHandle_p); |
1256 |
return; |
1257 |
} |
1258 |
isCellCentered[i_data]=FALSE; |
1259 |
break; |
1260 |
case FINLEY_NODES: |
1261 |
nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM; |
1262 |
if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) |
1263 |
{ |
1264 |
elementtype=FINLEY_ELEMENTS; |
1265 |
} |
1266 |
else |
1267 |
{ |
1268 |
Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file."); |
1269 |
fclose(fileHandle_p); |
1270 |
return; |
1271 |
} |
1272 |
isCellCentered[i_data]=FALSE; |
1273 |
break; |
1274 |
case FINLEY_ELEMENTS: |
1275 |
nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM; |
1276 |
if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) |
1277 |
{ |
1278 |
elementtype=FINLEY_ELEMENTS; |
1279 |
} |
1280 |
else |
1281 |
{ |
1282 |
Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file."); |
1283 |
fclose(fileHandle_p); |
1284 |
return; |
1285 |
} |
1286 |
isCellCentered[i_data]=TRUE; |
1287 |
break; |
1288 |
case FINLEY_FACE_ELEMENTS: |
1289 |
nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM; |
1290 |
if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_FACE_ELEMENTS) |
1291 |
{ |
1292 |
elementtype=FINLEY_FACE_ELEMENTS; |
1293 |
} |
1294 |
else |
1295 |
{ |
1296 |
Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file."); |
1297 |
fclose(fileHandle_p); |
1298 |
return; |
1299 |
} |
1300 |
isCellCentered[i_data]=TRUE; |
1301 |
break; |
1302 |
case FINLEY_POINTS: |
1303 |
nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM; |
1304 |
if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_POINTS) |
1305 |
{ |
1306 |
elementtype=FINLEY_POINTS; |
1307 |
} |
1308 |
else |
1309 |
{ |
1310 |
Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file."); |
1311 |
fclose(fileHandle_p); |
1312 |
return; |
1313 |
} |
1314 |
isCellCentered[i_data]=TRUE; |
1315 |
break; |
1316 |
case FINLEY_CONTACT_ELEMENTS_1: |
1317 |
nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM; |
1318 |
if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) |
1319 |
{ |
1320 |
elementtype=FINLEY_CONTACT_ELEMENTS_1; |
1321 |
} |
1322 |
else |
1323 |
{ |
1324 |
Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file."); |
1325 |
fclose(fileHandle_p); |
1326 |
return; |
1327 |
} |
1328 |
isCellCentered[i_data]=TRUE; |
1329 |
break; |
1330 |
case FINLEY_CONTACT_ELEMENTS_2: |
1331 |
nodetype = (nodetype == FINLEY_REDUCED_DEGREES_OF_FREEDOM) ? FINLEY_REDUCED_DEGREES_OF_FREEDOM : FINLEY_DEGREES_OF_FREEDOM; |
1332 |
if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) |
1333 |
{ |
1334 |
elementtype=FINLEY_CONTACT_ELEMENTS_1; |
1335 |
} |
1336 |
else |
1337 |
{ |
1338 |
Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file."); |
1339 |
fclose(fileHandle_p); |
1340 |
return; |
1341 |
} |
1342 |
isCellCentered[i_data]=TRUE; |
1343 |
break; |
1344 |
default: |
1345 |
sprintf(error_msg,"saveVTK: Finley does not know anything about function space type %d",getFunctionSpaceType(data_pp[i_data])); |
1346 |
Finley_setError(TYPE_ERROR,error_msg); |
1347 |
fclose(fileHandle_p); |
1348 |
return; |
1349 |
} |
1350 |
if (isCellCentered[i_data]) |
1351 |
{ |
1352 |
write_celldata=TRUE; |
1353 |
} |
1354 |
else |
1355 |
{ |
1356 |
write_pointdata=TRUE; |
1357 |
} |
1358 |
} |
1359 |
} |
1360 |
/* select nomber of points and the mesh component */ |
1361 |
numPoints = mesh_p->Nodes->numNodes; |
1362 |
if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) |
1363 |
{ |
1364 |
numPoints = mesh_p->Nodes->reducedNumNodes; |
1365 |
} |
1366 |
else |
1367 |
{ |
1368 |
numPoints = mesh_p->Nodes->numNodes; |
1369 |
} |
1370 |
if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS; |
1371 |
Finley_ElementFile* elements=NULL; |
1372 |
switch(elementtype) |
1373 |
{ |
1374 |
case FINLEY_ELEMENTS: |
1375 |
elements=mesh_p->Elements; |
1376 |
break; |
1377 |
case FINLEY_FACE_ELEMENTS: |
1378 |
elements=mesh_p->FaceElements; |
1379 |
break; |
1380 |
case FINLEY_POINTS: |
1381 |
elements=mesh_p->Points; |
1382 |
break; |
1383 |
case FINLEY_CONTACT_ELEMENTS_1: |
1384 |
elements=mesh_p->ContactElements; |
1385 |
break; |
1386 |
} |
1387 |
if (elements==NULL) |
1388 |
{ |
1389 |
Finley_setError(SYSTEM_ERROR,"saveVTK: undefined element file"); |
1390 |
fclose(fileHandle_p); |
1391 |
return; |
1392 |
} |
1393 |
/* map finley element type to VTK element type */ |
1394 |
numCells = elements->numElements; |
1395 |
ElementTypeId TypeId; |
1396 |
if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) |
1397 |
{ |
1398 |
TypeId = elements->LinearReferenceElement->Type->TypeId; |
1399 |
} |
1400 |
else |
1401 |
{ |
1402 |
TypeId = elements->ReferenceElement->Type->TypeId; |
1403 |
} |
1404 |
|
1405 |
switch(TypeId) |
1406 |
{ |
1407 |
case Point1: |
1408 |
case Line2Face: |
1409 |
case Line3Face: |
1410 |
case Point1_Contact: |
1411 |
case Line2Face_Contact: |
1412 |
case Line3Face_Contact: |
1413 |
cellType = VTK_VERTEX; |
1414 |
numVTKNodesPerElement = 1; |
1415 |
strcpy(elemTypeStr, "VTK_VERTEX"); |
1416 |
break; |
1417 |
|
1418 |
case Line2: |
1419 |
case Tri3Face: |
1420 |
case Rec4Face: |
1421 |
case Line2_Contact: |
1422 |
case Tri3_Contact: |
1423 |
case Tri3Face_Contact: |
1424 |
case Rec4Face_Contact: |
1425 |
cellType = VTK_LINE; |
1426 |
numVTKNodesPerElement = 2; |
1427 |
strcpy(elemTypeStr, "VTK_LINE"); |
1428 |
break; |
1429 |
|
1430 |
case Tri3: |
1431 |
case Tet4Face: |
1432 |
case Tet4Face_Contact: |
1433 |
cellType = VTK_TRIANGLE; |
1434 |
numVTKNodesPerElement = 3; |
1435 |
strcpy(elemTypeStr, "VTK_TRIANGLE"); |
1436 |
break; |
1437 |
|
1438 |
case Rec4: |
1439 |
case Hex8Face: |
1440 |
case Rec4_Contact: |
1441 |
case Hex8Face_Contact: |
1442 |
cellType = VTK_QUAD; |
1443 |
numVTKNodesPerElement = 4; |
1444 |
strcpy(elemTypeStr, "VTK_QUAD"); |
1445 |
break; |
1446 |
|
1447 |
case Tet4: |
1448 |
cellType = VTK_TETRA; |
1449 |
numVTKNodesPerElement = 4; |
1450 |
strcpy(elemTypeStr, "VTK_TETRA"); |
1451 |
break; |
1452 |
|
1453 |
case Hex8: |
1454 |
cellType = VTK_HEXAHEDRON; |
1455 |
numVTKNodesPerElement = 8; |
1456 |
strcpy(elemTypeStr, "VTK_HEXAHEDRON"); |
1457 |
break; |
1458 |
|
1459 |
case Line3: |
1460 |
case Tri6Face: |
1461 |
case Rec8Face: |
1462 |
case Line3_Contact: |
1463 |
case Tri6Face_Contact: |
1464 |
case Rec8Face_Contact: |
1465 |
cellType = VTK_QUADRATIC_EDGE; |
1466 |
numVTKNodesPerElement = 3; |
1467 |
strcpy(elemTypeStr, "VTK_QUADRATIC_EDGE"); |
1468 |
break; |
1469 |
|
1470 |
case Tri6: |
1471 |
case Tet10Face: |
1472 |
case Tri6_Contact: |
1473 |
case Tet10Face_Contact: |
1474 |
cellType = VTK_QUADRATIC_TRIANGLE; |
1475 |
numVTKNodesPerElement = 6; |
1476 |
strcpy(elemTypeStr, "VTK_QUADRATIC_TRIANGLE"); |
1477 |
break; |
1478 |
|
1479 |
case Rec8: |
1480 |
case Hex20Face: |
1481 |
case Rec8_Contact: |
1482 |
case Hex20Face_Contact: |
1483 |
cellType = VTK_QUADRATIC_QUAD; |
1484 |
numVTKNodesPerElement = 8; |
1485 |
strcpy(elemTypeStr, "VTK_QUADRATIC_QUAD"); |
1486 |
break; |
1487 |
|
1488 |
case Tet10: |
1489 |
cellType = VTK_QUADRATIC_TETRA; |
1490 |
numVTKNodesPerElement = 10; |
1491 |
strcpy(elemTypeStr, "VTK_QUADRATIC_TETRA"); |
1492 |
break; |
1493 |
|
1494 |
case Hex20: |
1495 |
cellType = VTK_QUADRATIC_HEXAHEDRON; |
1496 |
numVTKNodesPerElement = 20; |
1497 |
strcpy(elemTypeStr, "VTK_QUADRATIC_HEXAHEDRON"); |
1498 |
break; |
1499 |
|
1500 |
default: |
1501 |
sprintf(error_msg, "saveVTK: Element type %s is not supported by VTK",elements->ReferenceElement->Type->Name); |
1502 |
Finley_setError(VALUE_ERROR,error_msg); |
1503 |
fclose(fileHandle_p); |
1504 |
return; |
1505 |
} |
1506 |
/* xml header */ |
1507 |
fprintf(fileHandle_p, "<?xml version=\"1.0\"?>\n"); |
1508 |
fprintf(fileHandle_p, "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\">\n"); |
1509 |
|
1510 |
/* finley uses an unstructured mesh, so UnstructuredGrid *should* work */ |
1511 |
fprintf(fileHandle_p, "<UnstructuredGrid>\n"); |
1512 |
|
1513 |
/* is there only one "piece" to the data?? */ |
1514 |
fprintf(fileHandle_p, "<Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n",numPoints, numCells); |
1515 |
/* now for the points; equivalent to positions section in saveDX() */ |
1516 |
/* "The points element explicitly defines coordinates for each point |
1517 |
* individually. It contains one DataArray element describing an array |
1518 |
* with three components per value, each specifying the coordinates of one |
1519 |
* point" - from Vtk User's Guide |
1520 |
*/ |
1521 |
fprintf(fileHandle_p, "<Points>\n"); |
1522 |
/* |
1523 |
* the reason for this if statement is explained in the long comment below |
1524 |
*/ |
1525 |
nDim = mesh_p->Nodes->numDim; |
1526 |
fprintf(fileHandle_p, "<DataArray NumberOfComponents=\"%d\" type=\"Float32\" format=\"ascii\">\n",MAX(3,nDim)); |
1527 |
/* vtk/mayavi doesn't like 2D data, it likes 3D data with a degenerate |
1528 |
* third dimension to handle 2D data (like a sheet of paper). So, if |
1529 |
* nDim is 2, we have to append zeros to the array to get this third |
1530 |
* dimension, and keep the visualisers happy. |
1531 |
* Indeed, if nDim is less than 3, must pad all empty dimensions, so |
1532 |
* that the total number of dims is 3. |
1533 |
*/ |
1534 |
if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) |
1535 |
{ |
1536 |
for (i = 0; i < mesh_p->Nodes->numNodes; i++) |
1537 |
{ |
1538 |
if (mesh_p->Nodes->toReduced[i]>=0) |
1539 |
{ |
1540 |
for (j = 0; j < nDim; j++) fprintf(fileHandle_p, " %e",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]); |
1541 |
for (k=0; k<3-nDim; k++) fprintf(fileHandle_p, " %e",0.); |
1542 |
fprintf(fileHandle_p, "\n"); |
1543 |
} |
1544 |
} |
1545 |
} |
1546 |
else |
1547 |
{ |
1548 |
for (i = 0; i < mesh_p->Nodes->numNodes; i++) |
1549 |
{ |
1550 |
|
1551 |
for (j = 0; j < nDim; j++) fprintf(fileHandle_p, " %e",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]); |
1552 |
for (k=0; k<3-nDim; k++) fprintf(fileHandle_p, " %e",0.); |
1553 |
fprintf(fileHandle_p, "\n"); |
1554 |
} |
1555 |
|
1556 |
} |
1557 |
fprintf(fileHandle_p, "</DataArray>\n"); |
1558 |
fprintf(fileHandle_p, "</Points>\n"); |
1559 |
|
1560 |
/* write out the DataArray element for the connectivity */ |
1561 |
|
1562 |
int NN = elements->ReferenceElement->Type->numNodes; |
1563 |
fprintf(fileHandle_p, "<Cells>\n"); |
1564 |
fprintf(fileHandle_p, "<DataArray Name=\"connectivity\" type=\"Int32\" format=\"ascii\">\n"); |
1565 |
|
1566 |
if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) |
1567 |
{ |
1568 |
for (i = 0; i < numCells; i++) |
1569 |
{ |
1570 |
for (j = 0; j < numVTKNodesPerElement; j++) |
1571 |
fprintf(fileHandle_p,"%d ",mesh_p->Nodes->toReduced[elements->Nodes[INDEX2(elements->ReferenceElement->Type->linearNodes[j], i, NN)]]); |
1572 |
fprintf(fileHandle_p, "\n"); |
1573 |
} |
1574 |
} |
1575 |
else if (VTK_QUADRATIC_HEXAHEDRON==cellType) |
1576 |
{ |
1577 |
for (i = 0; i < numCells; i++) |
1578 |
{ |
1579 |
fprintf(fileHandle_p,"%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d\n", |
1580 |
elements->Nodes[INDEX2(0, i, NN)], |
1581 |
elements->Nodes[INDEX2(1, i, NN)], |
1582 |
elements->Nodes[INDEX2(2, i, NN)], |
1583 |
elements->Nodes[INDEX2(3, i, NN)], |
1584 |
elements->Nodes[INDEX2(4, i, NN)], |
1585 |
elements->Nodes[INDEX2(5, i, NN)], |
1586 |
elements->Nodes[INDEX2(6, i, NN)], |
1587 |
elements->Nodes[INDEX2(7, i, NN)], |
1588 |
elements->Nodes[INDEX2(8, i, NN)], |
1589 |
elements->Nodes[INDEX2(9, i, NN)], |
1590 |
elements->Nodes[INDEX2(10, i, NN)], |
1591 |
elements->Nodes[INDEX2(11, i, NN)], |
1592 |
elements->Nodes[INDEX2(16, i, NN)], |
1593 |
elements->Nodes[INDEX2(17, i, NN)], |
1594 |
elements->Nodes[INDEX2(18, i, NN)], |
1595 |
elements->Nodes[INDEX2(19, i, NN)], |
1596 |
elements->Nodes[INDEX2(12, i, NN)], |
1597 |
elements->Nodes[INDEX2(13, i, NN)], |
1598 |
elements->Nodes[INDEX2(14, i, NN)], |
1599 |
elements->Nodes[INDEX2(15, i, NN)]); |
1600 |
} |
1601 |
} |
1602 |
else if (numVTKNodesPerElement!=NN) |
1603 |
{ |
1604 |
for (i = 0; i < numCells; i++) |
1605 |
{ |
1606 |
for (j = 0; j < numVTKNodesPerElement; j++) fprintf(fileHandle_p,"%d ", elements->Nodes[INDEX2(elements->ReferenceElement->Type->geoNodes[j], i, NN)]); |
1607 |
fprintf(fileHandle_p, "\n"); |
1608 |
} |
1609 |
} |
1610 |
else |
1611 |
{ |
1612 |
for (i = 0; i < numCells; i++) |
1613 |
{ |
1614 |
for (j = 0; j < numVTKNodesPerElement; j++) fprintf(fileHandle_p,"%d ", elements->Nodes[INDEX2(j, i, NN)]); |
1615 |
fprintf(fileHandle_p, "\n"); |
1616 |
} |
1617 |
} |
1618 |
fprintf(fileHandle_p, "</DataArray>\n"); |
1619 |
|
1620 |
/* write out the DataArray element for the offsets */ |
1621 |
fprintf(fileHandle_p, "<DataArray Name=\"offsets\" type=\"Int32\" format=\"ascii\">\n"); |
1622 |
for (i=numVTKNodesPerElement; i<=numCells*numVTKNodesPerElement; i+=numVTKNodesPerElement) fprintf(fileHandle_p, "%d\n", i); |
1623 |
fprintf(fileHandle_p, "</DataArray>\n"); |
1624 |
|
1625 |
/* write out the DataArray element for the types */ |
1626 |
fprintf(fileHandle_p, "<DataArray Name=\"types\" type=\"UInt8\" format=\"ascii\">\n"); |
1627 |
for (i=0; i<numCells; i++) fprintf(fileHandle_p, "%d\n", cellType); |
1628 |
fprintf(fileHandle_p, "</DataArray>\n"); |
1629 |
|
1630 |
/* finish off the <Cells> element */ |
1631 |
fprintf(fileHandle_p, "</Cells>\n"); |
1632 |
|
1633 |
/* cell data */ |
1634 |
if (write_celldata) |
1635 |
{ |
1636 |
/* mark the active data arrays */ |
1637 |
bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE; |
1638 |
fprintf(fileHandle_p, "<CellData"); |
1639 |
for (i_data =0 ;i_data<num_data;++i_data) |
1640 |
{ |
1641 |
if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data]) |
1642 |
{ |
1643 |
/* if the rank == 0: --> scalar data |
1644 |
* if the rank == 1: --> vector data |
1645 |
* if the rank == 2: --> tensor data |
1646 |
*/ |
1647 |
switch(getDataPointRank(data_pp[i_data])) |
1648 |
{ |
1649 |
case 0: |
1650 |
if (! set_scalar) |
1651 |
{ |
1652 |
fprintf(fileHandle_p," Scalars=\"%s\"",names_p[i_data]); |
1653 |
set_scalar=TRUE; |
1654 |
} |
1655 |
break; |
1656 |
case 1: |
1657 |
if (! set_vector) |
1658 |
{ |
1659 |
fprintf(fileHandle_p," Vectors=\"%s\"",names_p[i_data]); |
1660 |
set_vector=TRUE; |
1661 |
} |
1662 |
break; |
1663 |
case 2: |
1664 |
if (! set_tensor) |
1665 |
{ |
1666 |
fprintf(fileHandle_p," Tensors=\"%s\"",names_p[i_data]); |
1667 |
set_tensor=TRUE; |
1668 |
} |
1669 |
break; |
1670 |
default: |
1671 |
sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]); |
1672 |
Finley_setError(VALUE_ERROR,error_msg); |
1673 |
fclose(fileHandle_p); |
1674 |
return; |
1675 |
} |
1676 |
} |
1677 |
} |
1678 |
fprintf(fileHandle_p, ">\n"); |
1679 |
/* write the arrays */ |
1680 |
for (i_data =0 ;i_data<num_data;++i_data) |
1681 |
{ |
1682 |
if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data]) |
1683 |
{ |
1684 |
numPointsPerSample = elements->ReferenceElement->numQuadNodes; |
1685 |
rank = getDataPointRank(data_pp[i_data]); |
1686 |
nComp = getDataPointSize(data_pp[i_data]); |
1687 |
nCompReqd=1; /* the number of components required by vtk */ |
1688 |
shape=0; |
1689 |
if (rank == 0) |
1690 |
{ |
1691 |
nCompReqd = 1; |
1692 |
} |
1693 |
else if (rank == 1) |
1694 |
{ |
1695 |
shape=getDataPointShape(data_pp[i_data], 0); |
1696 |
if (shape>3) |
1697 |
{ |
1698 |
Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components"); |
1699 |
fclose(fileHandle_p); |
1700 |
return; |
1701 |
} |
1702 |
nCompReqd = 3; |
1703 |
} |
1704 |
else |
1705 |
{ |
1706 |
shape=getDataPointShape(data_pp[i_data], 0); |
1707 |
if (shape>3 || shape != getDataPointShape(data_pp[i_data], 1)) |
1708 |
{ |
1709 |
Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape"); |
1710 |
fclose(fileHandle_p); |
1711 |
return; |
1712 |
} |
1713 |
nCompReqd = 9; |
1714 |
} |
1715 |
fprintf(fileHandle_p, "<DataArray Name=\"%s\" type=\"Float32\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd); |
1716 |
|
1717 |
double sampleAvg[nComp]; |
1718 |
for (i=0; i<numCells; i++) |
1719 |
{ |
1720 |
values = getSampleData(data_pp[i_data], i); |
1721 |
/* averaging over the number of points in the sample */ |
1722 |
for (k=0; k<nComp; k++) |
1723 |
{ |
1724 |
rtmp = 0.; |
1725 |
for (j=0; j<numPointsPerSample; j++) rtmp += values[INDEX2(k,j,nComp)]; |
1726 |
sampleAvg[k] = rtmp/numPointsPerSample; |
1727 |
} |
1728 |
/* if the number of required components is more than the number |
1729 |
* of actual components, pad with zeros |
1730 |
*/ |
1731 |
/* probably only need to get shape of first element */ |
1732 |
/* write the data different ways for scalar, vector and tensor */ |
1733 |
if (nCompReqd == 1) |
1734 |
{ |
1735 |
fprintf(fileHandle_p, " %e", sampleAvg[0]); |
1736 |
} |
1737 |
else if (nCompReqd == 3) |
1738 |
{ |
1739 |
/* write out the data */ |
1740 |
for (m=0; m<shape; m++) fprintf(fileHandle_p, " %e", sampleAvg[m]); |
1741 |
for (m=0; m<nCompReqd-shape; m++) fprintf(fileHandle_p, " %e", 0.); |
1742 |
} |
1743 |
else if (nCompReqd == 9) |
1744 |
{ |
1745 |
/* tensor data, so have a 3x3 matrix to output as a row |
1746 |
* of 9 data points */ |
1747 |
count = 0; |
1748 |
for (m=0; m<shape; m++) |
1749 |
{ |
1750 |
for (n=0; n<shape; n++) |
1751 |
{ |
1752 |
fprintf(fileHandle_p, " %e", sampleAvg[count]); |
1753 |
count++; |
1754 |
} |
1755 |
for (n=0; n<3-shape; n++) fprintf(fileHandle_p, " %e", 0.); |
1756 |
} |
1757 |
for (m=0; m<3-shape; m++) |
1758 |
for (n=0; n<3; n++) fprintf(fileHandle_p, " %e", 0.); |
1759 |
} |
1760 |
fprintf(fileHandle_p, "\n"); |
1761 |
} |
1762 |
fprintf(fileHandle_p, "</DataArray>\n"); |
1763 |
} |
1764 |
} |
1765 |
fprintf(fileHandle_p, "</CellData>\n"); |
1766 |
} |
1767 |
/* point data */ |
1768 |
if (write_pointdata) |
1769 |
{ |
1770 |
/* mark the active data arrays */ |
1771 |
bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE; |
1772 |
fprintf(fileHandle_p, "<PointData"); |
1773 |
for (i_data =0 ;i_data<num_data;++i_data) |
1774 |
{ |
1775 |
if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data]) |
1776 |
{ |
1777 |
/* if the rank == 0: --> scalar data |
1778 |
* if the rank == 1: --> vector data |
1779 |
* if the rank == 2: --> tensor data |
1780 |
*/ |
1781 |
switch(getDataPointRank(data_pp[i_data])) |
1782 |
{ |
1783 |
case 0: |
1784 |
if (! set_scalar) |
1785 |
{ |
1786 |
fprintf(fileHandle_p," Scalars=\"%s\"",names_p[i_data]); |
1787 |
set_scalar=TRUE; |
1788 |
} |
1789 |
break; |
1790 |
case 1: |
1791 |
if (! set_vector) |
1792 |
{ |
1793 |
fprintf(fileHandle_p," Vectors=\"%s\"",names_p[i_data]); |
1794 |
set_vector=TRUE; |
1795 |
} |
1796 |
break; |
1797 |
case 2: |
1798 |
if (! set_tensor) |
1799 |
{ |
1800 |
fprintf(fileHandle_p," Tensors=\"%s\"",names_p[i_data]); |
1801 |
set_tensor=TRUE; |
1802 |
} |
1803 |
break; |
1804 |
default: |
1805 |
sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]); |
1806 |
Finley_setError(VALUE_ERROR,error_msg); |
1807 |
fclose(fileHandle_p); |
1808 |
return; |
1809 |
} |
1810 |
} |
1811 |
} |
1812 |
fprintf(fileHandle_p, ">\n"); |
1813 |
/* write the arrays */ |
1814 |
for (i_data =0 ;i_data<num_data;++i_data) |
1815 |
{ |
1816 |
if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data]) |
1817 |
{ |
1818 |
numPointsPerSample = elements->ReferenceElement->numQuadNodes; |
1819 |
rank = getDataPointRank(data_pp[i_data]); |
1820 |
nComp = getDataPointSize(data_pp[i_data]); |
1821 |
nCompReqd=1; /* the number of components required by vtk */ |
1822 |
shape=0; |
1823 |
if (rank == 0) |
1824 |
{ |
1825 |
nCompReqd = 1; |
1826 |
} |
1827 |
else if (rank == 1) |
1828 |
{ |
1829 |
shape=getDataPointShape(data_pp[i_data], 0); |
1830 |
if (shape>3) |
1831 |
{ |
1832 |
Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components"); |
1833 |
fclose(fileHandle_p); |
1834 |
return; |
1835 |
} |
1836 |
nCompReqd = 3; |
1837 |
} |
1838 |
else |
1839 |
{ |
1840 |
shape=getDataPointShape(data_pp[i_data], 0); |
1841 |
if (shape>3 || shape != getDataPointShape(data_pp[i_data], 1)) |
1842 |
{ |
1843 |
Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape"); |
1844 |
fclose(fileHandle_p); |
1845 |
return; |
1846 |
} |
1847 |
nCompReqd = 9; |
1848 |
} |
1849 |
fprintf(fileHandle_p, "<DataArray Name=\"%s\" type=\"Float32\" NumberOfComponents=\"%d\" format=\"ascii\">\n",names_p[i_data], nCompReqd); |
1850 |
/* write out the data */ |
1851 |
/* if the number of required components is more than the number |
1852 |
* of actual components, pad with zeros |
1853 |
*/ |
1854 |
bool_t do_write=TRUE; |
1855 |
for (i=0; i<mesh_p->Nodes->numNodes; i++) |
1856 |
{ |
1857 |
if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) |
1858 |
{ |
1859 |
if (mesh_p->Nodes->toReduced[i]>=0) |
1860 |
{ |
1861 |
switch(getFunctionSpaceType(data_pp[i_data])) |
1862 |
{ |
1863 |
case FINLEY_DEGREES_OF_FREEDOM: |
1864 |
values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]); |
1865 |
break; |
1866 |
case FINLEY_REDUCED_DEGREES_OF_FREEDOM: |
1867 |
values = getSampleData(data_pp[i_data],mesh_p->Nodes->reducedDegreeOfFreedom[i]); |
1868 |
break; |
1869 |
case FINLEY_NODES: |
1870 |
values = getSampleData(data_pp[i_data],i); |
1871 |
break; |
1872 |
} |
1873 |
do_write=TRUE; |
1874 |
} |
1875 |
else |
1876 |
{ |
1877 |
do_write=FALSE; |
1878 |
} |
1879 |
} |
1880 |
else |
1881 |
{ |
1882 |
do_write=TRUE; |
1883 |
switch(getFunctionSpaceType(data_pp[i_data])) |
1884 |
{ |
1885 |
case FINLEY_DEGREES_OF_FREEDOM: |
1886 |
values = getSampleData(data_pp[i_data],mesh_p->Nodes->degreeOfFreedom[i]); |
1887 |
break; |
1888 |
case FINLEY_NODES: |
1889 |
values = getSampleData(data_pp[i_data],i); |
1890 |
break; |
1891 |
} |
1892 |
} |
1893 |
/* write the data different ways for scalar, vector and tensor */ |
1894 |
if (do_write) |
1895 |
{ |
1896 |
if (nCompReqd == 1) |
1897 |
{ |
1898 |
fprintf(fileHandle_p, " %e", values[0]); |
1899 |
} |
1900 |
else if (nCompReqd == 3) |
1901 |
{ |
1902 |
for (m=0; m<shape; m++) fprintf(fileHandle_p, " %e", values[m]); |
1903 |
for (m=0; m<nCompReqd-shape; m++) fprintf(fileHandle_p, " %e", 0.); |
1904 |
} |
1905 |
else if (nCompReqd == 9) |
1906 |
{ |
1907 |
/* tensor data, so have a 3x3 matrix to output as a row |
1908 |
* of 9 data points */ |
1909 |
count = 0; |
1910 |
for (m=0; m<shape; m++) |
1911 |
{ |
1912 |
for (n=0; n<shape; n++) |
1913 |
{ |
1914 |
fprintf(fileHandle_p, " %e", values[count]); |
1915 |
count++; |
1916 |
} |
1917 |
for (n=0; n<3-shape; n++) fprintf(fileHandle_p, " %e", 0.); |
1918 |
} |
1919 |
for (m=0; m<3-shape; m++) |
1920 |
for (n=0; n<3; n++) fprintf(fileHandle_p, " %e", 0.); |
1921 |
} |
1922 |
fprintf(fileHandle_p, "\n"); |
1923 |
} |
1924 |
} |
1925 |
fprintf(fileHandle_p, "</DataArray>\n"); |
1926 |
} |
1927 |
} |
1928 |
fprintf(fileHandle_p, "</PointData>\n"); |
1929 |
} |
1930 |
/* finish off the piece */ |
1931 |
fprintf(fileHandle_p, "</Piece>\n"); |
1932 |
|
1933 |
fprintf(fileHandle_p, "</UnstructuredGrid>\n"); |
1934 |
/* write the xml footer */ |
1935 |
fprintf(fileHandle_p, "</VTKFile>\n"); |
1936 |
/* close the file */ |
1937 |
fclose(fileHandle_p); |
1938 |
return; |
1939 |
} |
1940 |
#endif |
1941 |
|