Actual source code: plexvtu.c
petsc-3.4.5 2014-06-29
1: #include <petsc-private/dmpleximpl.h>
2: #include <../src/sys/classes/viewer/impls/vtk/vtkvimpl.h>
4: typedef struct {
5: PetscInt nvertices;
6: PetscInt ncells;
7: PetscInt nconn; /* number of entries in cell->vertex connectivity array */
8: } PieceInfo;
10: #if defined(PETSC_USE_REAL_SINGLE)
11: static const char precision[] = "Float32";
12: #elif defined(PETSC_USE_REAL_DOUBLE)
13: static const char precision[] = "Float64";
14: #else
15: static const char precision[] = "UnknownPrecision";
16: #endif
20: static PetscErrorCode TransferWrite(PetscViewer viewer,FILE *fp,PetscMPIInt srank,PetscMPIInt root,const void *send,void *recv,PetscMPIInt count,PetscDataType datatype,PetscMPIInt tag)
21: {
22: PetscMPIInt rank;
24: MPI_Comm comm;
25: MPI_Datatype mpidatatype;
28: PetscObjectGetComm((PetscObject)viewer,&comm);
29: MPI_Comm_rank(comm,&rank);
30: PetscDataTypeToMPIDataType(datatype,&mpidatatype);
32: if (rank == srank && rank != root) {
33: MPI_Send((void*)send,count,mpidatatype,root,tag,comm);
34: } else if (rank == root) {
35: const void *buffer;
36: if (root == srank) { /* self */
37: buffer = send;
38: } else {
39: MPI_Status status;
40: PetscMPIInt nrecv;
41: MPI_Recv(recv,count,mpidatatype,srank,tag,comm,&status);
42: MPI_Get_count(&status,mpidatatype,&nrecv);
43: if (count != nrecv) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Array size mismatch");
44: buffer = recv;
45: }
46: PetscViewerVTKFWrite(viewer,fp,buffer,count,datatype);
47: }
48: return(0);
49: }
53: static PetscErrorCode DMPlexGetVTKConnectivity(DM dm,PieceInfo *piece,PetscVTKInt **oconn,PetscVTKInt **ooffsets,PetscVTKType **otypes)
54: {
56: PetscVTKInt *conn,*offsets;
57: PetscVTKType *types;
58: PetscInt dim,vStart,vEnd,cStart,cEnd,pStart,pEnd,cellHeight,cMax,numLabelCells,hasLabel,c,v,countcell,countconn;
61: PetscMalloc3(piece->nconn,PetscVTKInt,&conn,piece->ncells,PetscVTKInt,&offsets,piece->ncells,PetscVTKType,&types);
63: DMPlexGetDimension(dm,&dim);
64: DMPlexGetChart(dm,&pStart,&pEnd);
65: DMPlexGetVTKCellHeight(dm, &cellHeight);
66: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
67: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
68: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
69: if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
70: DMPlexGetStratumSize(dm, "vtk", 1, &numLabelCells);
71: hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;
73: countcell = 0;
74: countconn = 0;
75: for (c = cStart; c < cEnd; ++c) {
76: PetscInt *closure = NULL;
77: PetscInt closureSize,nverts,celltype,startoffset;
79: if (hasLabel) {
80: PetscInt value;
82: DMPlexGetLabelValue(dm, "vtk", c, &value);
83: if (value != 1) continue;
84: }
85: startoffset = countconn;
86: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
87: for (v = 0; v < closureSize*2; v += 2) {
88: if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
89: conn[countconn++] = closure[v] - vStart;
90: }
91: }
92: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
94: offsets[countcell] = countconn;
96: nverts = countconn - startoffset;
97: DMPlexVTKGetCellType(dm,dim,nverts,&celltype);
99: types[countcell] = celltype;
100: countcell++;
101: }
102: if (countcell != piece->ncells) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent cell count");
103: if (countconn != piece->nconn) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent connectivity count");
104: *oconn = conn;
105: *ooffsets = offsets;
106: *otypes = types;
107: return(0);
108: }
112: /*
113: Write all fields that have been provided to the viewer
114: Multi-block XML format with binary appended data.
115: */
116: PetscErrorCode DMPlexVTKWriteAll_VTU(DM dm,PetscViewer viewer)
117: {
118: MPI_Comm comm;
119: PetscViewer_VTK *vtk = (PetscViewer_VTK*)viewer->data;
120: PetscViewerVTKObjectLink link;
121: FILE *fp;
122: PetscMPIInt rank,size,tag;
123: PetscErrorCode ierr;
124: PetscInt dim,cellHeight,cStart,cEnd,vStart,vEnd,cMax,numLabelCells,hasLabel,c,v,r,i;
125: PieceInfo piece,*gpiece = NULL;
126: void *buffer = NULL;
129: PetscObjectGetComm((PetscObject)dm,&comm);
130: #if defined(PETSC_USE_COMPLEX)
131: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Complex values not supported");
132: #endif
133: MPI_Comm_size(comm,&size);
134: MPI_Comm_rank(comm,&rank);
135: tag = ((PetscObject)viewer)->tag;
137: PetscFOpen(comm,vtk->filename,"wb",&fp);
138: PetscFPrintf(comm,fp,"<?xml version=\"1.0\"?>\n");
139: #if defined(PETSC_WORDS_BIGENDIAN)
140: PetscFPrintf(comm,fp,"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");
141: #else
142: PetscFPrintf(comm,fp,"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");
143: #endif
144: PetscFPrintf(comm,fp," <UnstructuredGrid>\n");
146: DMPlexGetDimension(dm, &dim);
147: DMPlexGetVTKCellHeight(dm, &cellHeight);
148: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
149: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
150: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
151: if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
152: DMPlexGetStratumSize(dm, "vtk", 1, &numLabelCells);
154: hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;
155: piece.nvertices = vEnd - vStart;
156: piece.ncells = 0;
157: piece.nconn = 0;
158: for (c = cStart; c < cEnd; ++c) {
159: PetscInt *closure = NULL;
160: PetscInt closureSize;
162: if (hasLabel) {
163: PetscInt value;
165: DMPlexGetLabelValue(dm, "vtk", c, &value);
166: if (value != 1) continue;
167: }
168: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
169: for (v = 0; v < closureSize*2; v += 2) {
170: if ((closure[v] >= vStart) && (closure[v] < vEnd)) piece.nconn++;
171: }
172: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
173: piece.ncells++;
174: }
175: if (!rank) {PetscMalloc(size*sizeof(piece),&gpiece);}
176: MPI_Gather((PetscInt*)&piece,sizeof(piece)/sizeof(PetscInt),MPIU_INT,(PetscInt*)gpiece,sizeof(piece)/sizeof(PetscInt),MPIU_INT,0,comm);
178: /*
179: * Write file header
180: */
181: if (!rank) {
182: PetscInt boffset = 0;
184: for (r=0; r<size; r++) {
185: PetscFPrintf(PETSC_COMM_SELF,fp," <Piece NumberOfPoints=\"%D\" NumberOfCells=\"%D\">\n",gpiece[r].nvertices,gpiece[r].ncells);
186: /* Coordinate positions */
187: PetscFPrintf(PETSC_COMM_SELF,fp," <Points>\n");
188: PetscFPrintf(PETSC_COMM_SELF,fp," <DataArray type=\"%s\" Name=\"Position\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,boffset);
189: boffset += gpiece[r].nvertices*3*sizeof(PetscScalar) + sizeof(int);
190: PetscFPrintf(PETSC_COMM_SELF,fp," </Points>\n");
191: /* Cell connectivity */
192: PetscFPrintf(PETSC_COMM_SELF,fp," <Cells>\n");
193: PetscFPrintf(PETSC_COMM_SELF,fp," <DataArray type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",boffset);
194: boffset += gpiece[r].nconn*sizeof(PetscInt) + sizeof(int);
195: PetscFPrintf(PETSC_COMM_SELF,fp," <DataArray type=\"Int32\" Name=\"offsets\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",boffset);
196: boffset += gpiece[r].ncells*sizeof(PetscInt) + sizeof(int);
197: PetscFPrintf(PETSC_COMM_SELF,fp," <DataArray type=\"UInt8\" Name=\"types\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",boffset);
198: boffset += gpiece[r].ncells*sizeof(unsigned char) + sizeof(int);
199: PetscFPrintf(PETSC_COMM_SELF,fp," </Cells>\n");
201: /*
202: * Cell Data headers
203: */
204: PetscFPrintf(PETSC_COMM_SELF,fp," <CellData>\n");
205: PetscFPrintf(PETSC_COMM_SELF,fp," <DataArray type=\"Int32\" Name=\"Rank\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",boffset);
206: boffset += gpiece[r].ncells*sizeof(int) + sizeof(int);
207: /* all the vectors */
208: for (link=vtk->link; link; link=link->next) {
209: Vec X = (Vec)link->vec;
210: PetscInt bs,nfields,field;
211: const char *vecname = "";
212: if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
213: if (((PetscObject)X)->name || link != vtk->link) { /* If the object is already named, use it. If it is past the first link, name it to disambiguate. */
214: PetscObjectGetName((PetscObject)X,&vecname);
215: }
216: PetscSectionGetDof(dm->defaultSection,cStart,&bs);
217: PetscSectionGetNumFields(dm->defaultSection,&nfields);
218: for (field=0,i=0; field<(nfields?nfields:1); field++) {
219: PetscInt fbs,j;
220: const char *fieldname = NULL;
221: char buf[256];
222: if (nfields) { /* We have user-defined fields/components */
223: PetscSectionGetFieldDof(dm->defaultSection,cStart,field,&fbs);
224: PetscSectionGetFieldName(dm->defaultSection,field,&fieldname);
225: } else fbs = bs; /* Say we have one field with 'bs' components */
226: if (!fieldname) {
227: PetscSNPrintf(buf,sizeof(buf),"CellField%D",field);
228: fieldname = buf;
229: }
230: for (j=0; j<fbs; j++) {
231: PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s.%D\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,j,boffset);
232: boffset += gpiece[r].ncells*sizeof(PetscScalar) + sizeof(int);
233: i++;
234: }
235: }
236: if (i != bs) SETERRQ2(comm,PETSC_ERR_PLIB,"Total number of field components %D != block size %D",i,bs);
237: }
238: PetscFPrintf(PETSC_COMM_SELF,fp," </CellData>\n");
240: /*
241: * Point Data headers
242: */
243: PetscFPrintf(PETSC_COMM_SELF,fp," <PointData>\n");
244: for (link=vtk->link; link; link=link->next) {
245: Vec X = (Vec)link->vec;
246: PetscInt bs,nfields,field;
247: const char *vecname = "";
248: if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
249: if (((PetscObject)X)->name || link != vtk->link) { /* If the object is already named, use it. If it is past the first link, name it to disambiguate. */
250: PetscObjectGetName((PetscObject)X,&vecname);
251: }
252: PetscSectionGetDof(dm->defaultSection,vStart,&bs);
253: PetscSectionGetNumFields(dm->defaultSection,&nfields);
254: for (field=0,i=0; field<(nfields?nfields:1); field++) {
255: PetscInt fbs,j;
256: const char *fieldname = NULL;
257: char buf[256];
258: if (nfields) { /* We have user-defined fields/components */
259: PetscSectionGetFieldDof(dm->defaultSection,vStart,field,&fbs);
260: PetscSectionGetFieldName(dm->defaultSection,field,&fieldname);
261: } else fbs = bs; /* Say we have one field with 'bs' components */
262: if (!fieldname) {
263: PetscSNPrintf(buf,sizeof(buf),"PointField%D",field);
264: fieldname = buf;
265: }
266: for (j=0; j<fbs; j++) {
267: PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s.%D\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,j,boffset);
268: boffset += gpiece[r].nvertices*sizeof(PetscScalar) + sizeof(int);
269: }
270: }
271: }
272: PetscFPrintf(PETSC_COMM_SELF,fp," </PointData>\n");
274: PetscFPrintf(PETSC_COMM_SELF,fp," </Piece>\n");
275: }
276: }
278: PetscFPrintf(comm,fp," </UnstructuredGrid>\n");
279: PetscFPrintf(comm,fp," <AppendedData encoding=\"raw\">\n");
280: PetscFPrintf(comm,fp,"_");
282: if (!rank) {
283: PetscInt maxsize = 0;
284: for (r=0; r<size; r++) {
285: maxsize = PetscMax(maxsize, (PetscInt) (gpiece[r].nvertices*3*sizeof(PetscScalar)));
286: maxsize = PetscMax(maxsize, (PetscInt) (gpiece[r].ncells*sizeof(PetscScalar)));
287: maxsize = PetscMax(maxsize, (PetscInt) (gpiece[r].nconn*sizeof(PetscVTKInt)));
288: }
289: PetscMalloc(maxsize,&buffer);
290: }
291: for (r=0; r<size; r++) {
292: if (r == rank) {
293: PetscInt nsend;
294: { /* Position */
295: const PetscScalar *x;
296: PetscScalar *y = NULL;
297: Vec coords;
298: nsend = piece.nvertices*3;
299: DMGetCoordinatesLocal(dm,&coords);
300: VecGetArrayRead(coords,&x);
301: if (dim != 3) {
302: PetscMalloc(piece.nvertices*3*sizeof(PetscScalar),&y);
303: for (i=0; i<piece.nvertices; i++) {
304: y[i*3+0] = x[i*dim+0];
305: y[i*3+1] = (dim > 1) ? x[i*dim+1] : 0;
306: y[i*3+2] = 0;
307: }
308: }
309: TransferWrite(viewer,fp,r,0,y ? y : x,buffer,nsend,PETSC_SCALAR,tag);
310: PetscFree(y);
311: VecRestoreArrayRead(coords,&x);
312: }
313: { /* Connectivity, offsets, types */
314: PetscVTKInt *connectivity = NULL,*offsets;
315: PetscVTKType *types;
316: DMPlexGetVTKConnectivity(dm,&piece,&connectivity,&offsets,&types);
317: TransferWrite(viewer,fp,r,0,connectivity,buffer,piece.nconn,PETSC_INT32,tag);
318: TransferWrite(viewer,fp,r,0,offsets,buffer,piece.ncells,PETSC_INT32,tag);
319: TransferWrite(viewer,fp,r,0,types,buffer,piece.ncells,PETSC_UINT8,tag);
320: PetscFree3(connectivity,offsets,types);
321: }
322: { /* Owners (cell data) */
323: PetscVTKInt *owners;
324: PetscMalloc(piece.ncells*sizeof(PetscVTKInt),&owners);
325: for (i=0; i<piece.ncells; i++) owners[i] = rank;
326: TransferWrite(viewer,fp,r,0,owners,buffer,piece.ncells,PETSC_INT32,tag);
327: PetscFree(owners);
328: }
329: /* Cell data */
330: for (link=vtk->link; link; link=link->next) {
331: Vec X = (Vec)link->vec;
332: const PetscScalar *x;
333: PetscScalar *y;
334: PetscInt bs;
335: if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
336: PetscSectionGetDof(dm->defaultSection,cStart,&bs);
337: VecGetArrayRead(X,&x);
338: PetscMalloc(piece.ncells*sizeof(PetscScalar),&y);
339: for (i=0; i<bs; i++) {
340: PetscInt cnt;
341: for (c=cStart,cnt=0; c<cEnd; c++) {
342: const PetscScalar *xpoint;
343: if (hasLabel) { /* Ignore some cells */
344: PetscInt value;
345: DMPlexGetLabelValue(dm, "vtk", c, &value);
346: if (value != 1) continue;
347: }
348: DMPlexPointLocalRead(dm,c,x,&xpoint);
349: y[cnt++] = xpoint[i];
350: }
351: if (cnt != piece.ncells) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not match");
352: TransferWrite(viewer,fp,r,0,y,buffer,piece.ncells,PETSC_SCALAR,tag);
353: }
354: PetscFree(y);
355: VecRestoreArrayRead(X,&x);
356: }
358: for (link=vtk->link; link; link=link->next) {
359: Vec X = (Vec)link->vec;
360: const PetscScalar *x;
361: PetscScalar *y;
362: PetscInt bs;
363: if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
364: PetscSectionGetDof(dm->defaultSection,vStart,&bs);
365: VecGetArrayRead(X,&x);
366: PetscMalloc(piece.nvertices*sizeof(PetscScalar),&y);
367: for (i=0; i<bs; i++) {
368: PetscInt cnt;
369: for (v=vStart,cnt=0; v<vEnd; v++) {
370: const PetscScalar *xpoint;
371: DMPlexPointLocalRead(dm,v,x,&xpoint);
372: y[cnt++] = xpoint[i];
373: }
374: if (cnt != piece.nvertices) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not match");
375: TransferWrite(viewer,fp,r,0,y,buffer,piece.nvertices,PETSC_SCALAR,tag);
376: }
377: PetscFree(y);
378: VecRestoreArrayRead(X,&x);
379: }
380: } else if (!rank) {
381: TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].nvertices*3,PETSC_SCALAR,tag); /* positions */
382: TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].nconn,PETSC_INT32,tag); /* connectivity */
383: TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,PETSC_INT32,tag); /* offsets */
384: TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,PETSC_UINT8,tag); /* types */
385: TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,PETSC_INT32,tag); /* owner rank (cells) */
386: /* all cell data */
387: for (link=vtk->link; link; link=link->next) {
388: PetscInt bs;
389: if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
390: PetscSectionGetDof(dm->defaultSection,cStart,&bs);
391: for (i=0; i<bs; i++) {
392: TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,PETSC_SCALAR,tag);
393: }
394: }
395: /* all point data */
396: for (link=vtk->link; link; link=link->next) {
397: PetscInt bs;
398: if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
399: PetscSectionGetDof(dm->defaultSection,vStart,&bs);
400: for (i=0; i<bs; i++) {
401: TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].nvertices,PETSC_SCALAR,tag);
402: }
403: }
404: }
405: }
406: PetscFree(gpiece);
407: PetscFree(buffer);
408: PetscFPrintf(comm,fp,"\n </AppendedData>\n");
409: PetscFPrintf(comm,fp,"</VTKFile>\n");
410: return(0);
411: }