Actual source code: plexvtu.c
petsc-3.5.4 2015-05-23
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,&conn,piece->ncells,&offsets,piece->ncells,&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,nC=0;
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: ++nC;
91: }
92: }
93: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
94: DMPlexInvertCell(dim, nC, &conn[countconn-nC]);
96: offsets[countcell] = countconn;
98: nverts = countconn - startoffset;
99: DMPlexVTKGetCellType(dm,dim,nverts,&celltype);
101: types[countcell] = celltype;
102: countcell++;
103: }
104: if (countcell != piece->ncells) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent cell count");
105: if (countconn != piece->nconn) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent connectivity count");
106: *oconn = conn;
107: *ooffsets = offsets;
108: *otypes = types;
109: return(0);
110: }
114: /*
115: Write all fields that have been provided to the viewer
116: Multi-block XML format with binary appended data.
117: */
118: PetscErrorCode DMPlexVTKWriteAll_VTU(DM dm,PetscViewer viewer)
119: {
120: MPI_Comm comm;
121: PetscViewer_VTK *vtk = (PetscViewer_VTK*)viewer->data;
122: PetscViewerVTKObjectLink link;
123: FILE *fp;
124: PetscMPIInt rank,size,tag;
125: PetscErrorCode ierr;
126: PetscInt dim,cellHeight,cStart,cEnd,vStart,vEnd,cMax,numLabelCells,hasLabel,c,v,r,i;
127: PieceInfo piece,*gpiece = NULL;
128: void *buffer = NULL;
131: PetscObjectGetComm((PetscObject)dm,&comm);
132: #if defined(PETSC_USE_COMPLEX)
133: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Complex values not supported");
134: #endif
135: MPI_Comm_size(comm,&size);
136: MPI_Comm_rank(comm,&rank);
137: tag = ((PetscObject)viewer)->tag;
139: PetscFOpen(comm,vtk->filename,"wb",&fp);
140: PetscFPrintf(comm,fp,"<?xml version=\"1.0\"?>\n");
141: #if defined(PETSC_WORDS_BIGENDIAN)
142: PetscFPrintf(comm,fp,"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");
143: #else
144: PetscFPrintf(comm,fp,"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");
145: #endif
146: PetscFPrintf(comm,fp," <UnstructuredGrid>\n");
148: DMPlexGetDimension(dm, &dim);
149: DMPlexGetVTKCellHeight(dm, &cellHeight);
150: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
151: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
152: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
153: if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
154: DMPlexGetStratumSize(dm, "vtk", 1, &numLabelCells);
156: hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;
157: piece.nvertices = vEnd - vStart;
158: piece.ncells = 0;
159: piece.nconn = 0;
160: for (c = cStart; c < cEnd; ++c) {
161: PetscInt *closure = NULL;
162: PetscInt closureSize;
164: if (hasLabel) {
165: PetscInt value;
167: DMPlexGetLabelValue(dm, "vtk", c, &value);
168: if (value != 1) continue;
169: }
170: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
171: for (v = 0; v < closureSize*2; v += 2) {
172: if ((closure[v] >= vStart) && (closure[v] < vEnd)) piece.nconn++;
173: }
174: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
175: piece.ncells++;
176: }
177: if (!rank) {PetscMalloc1(size,&gpiece);}
178: MPI_Gather((PetscInt*)&piece,sizeof(piece)/sizeof(PetscInt),MPIU_INT,(PetscInt*)gpiece,sizeof(piece)/sizeof(PetscInt),MPIU_INT,0,comm);
180: /*
181: * Write file header
182: */
183: if (!rank) {
184: PetscInt boffset = 0;
186: for (r=0; r<size; r++) {
187: PetscFPrintf(PETSC_COMM_SELF,fp," <Piece NumberOfPoints=\"%D\" NumberOfCells=\"%D\">\n",gpiece[r].nvertices,gpiece[r].ncells);
188: /* Coordinate positions */
189: PetscFPrintf(PETSC_COMM_SELF,fp," <Points>\n");
190: PetscFPrintf(PETSC_COMM_SELF,fp," <DataArray type=\"%s\" Name=\"Position\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,boffset);
191: boffset += gpiece[r].nvertices*3*sizeof(PetscScalar) + sizeof(int);
192: PetscFPrintf(PETSC_COMM_SELF,fp," </Points>\n");
193: /* Cell connectivity */
194: PetscFPrintf(PETSC_COMM_SELF,fp," <Cells>\n");
195: PetscFPrintf(PETSC_COMM_SELF,fp," <DataArray type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",boffset);
196: boffset += gpiece[r].nconn*sizeof(PetscInt) + sizeof(int);
197: PetscFPrintf(PETSC_COMM_SELF,fp," <DataArray type=\"Int32\" Name=\"offsets\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",boffset);
198: boffset += gpiece[r].ncells*sizeof(PetscInt) + sizeof(int);
199: PetscFPrintf(PETSC_COMM_SELF,fp," <DataArray type=\"UInt8\" Name=\"types\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",boffset);
200: boffset += gpiece[r].ncells*sizeof(unsigned char) + sizeof(int);
201: PetscFPrintf(PETSC_COMM_SELF,fp," </Cells>\n");
203: /*
204: * Cell Data headers
205: */
206: PetscFPrintf(PETSC_COMM_SELF,fp," <CellData>\n");
207: PetscFPrintf(PETSC_COMM_SELF,fp," <DataArray type=\"Int32\" Name=\"Rank\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",boffset);
208: boffset += gpiece[r].ncells*sizeof(int) + sizeof(int);
209: /* all the vectors */
210: for (link=vtk->link; link; link=link->next) {
211: Vec X = (Vec)link->vec;
212: PetscInt bs,nfields,field;
213: const char *vecname = "";
214: if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
215: 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. */
216: PetscObjectGetName((PetscObject)X,&vecname);
217: }
218: PetscSectionGetDof(dm->defaultSection,cStart,&bs);
219: PetscSectionGetNumFields(dm->defaultSection,&nfields);
220: for (field=0,i=0; field<(nfields?nfields:1); field++) {
221: PetscInt fbs,j;
222: const char *fieldname = NULL;
223: char buf[256];
224: if (nfields) { /* We have user-defined fields/components */
225: PetscSectionGetFieldDof(dm->defaultSection,cStart,field,&fbs);
226: PetscSectionGetFieldName(dm->defaultSection,field,&fieldname);
227: } else fbs = bs; /* Say we have one field with 'bs' components */
228: if (!fieldname) {
229: PetscSNPrintf(buf,sizeof(buf),"CellField%D",field);
230: fieldname = buf;
231: }
232: for (j=0; j<fbs; j++) {
233: PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s.%D\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,j,boffset);
234: boffset += gpiece[r].ncells*sizeof(PetscScalar) + sizeof(int);
235: i++;
236: }
237: }
238: if (i != bs) SETERRQ2(comm,PETSC_ERR_PLIB,"Total number of field components %D != block size %D",i,bs);
239: }
240: PetscFPrintf(PETSC_COMM_SELF,fp," </CellData>\n");
242: /*
243: * Point Data headers
244: */
245: PetscFPrintf(PETSC_COMM_SELF,fp," <PointData>\n");
246: for (link=vtk->link; link; link=link->next) {
247: Vec X = (Vec)link->vec;
248: PetscInt bs,nfields,field;
249: const char *vecname = "";
250: if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
251: 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. */
252: PetscObjectGetName((PetscObject)X,&vecname);
253: }
254: PetscSectionGetDof(dm->defaultSection,vStart,&bs);
255: PetscSectionGetNumFields(dm->defaultSection,&nfields);
256: for (field=0,i=0; field<(nfields?nfields:1); field++) {
257: PetscInt fbs,j;
258: const char *fieldname = NULL;
259: char buf[256];
260: if (nfields) { /* We have user-defined fields/components */
261: PetscSectionGetFieldDof(dm->defaultSection,vStart,field,&fbs);
262: PetscSectionGetFieldName(dm->defaultSection,field,&fieldname);
263: } else fbs = bs; /* Say we have one field with 'bs' components */
264: if (!fieldname) {
265: PetscSNPrintf(buf,sizeof(buf),"PointField%D",field);
266: fieldname = buf;
267: }
268: for (j=0; j<fbs; j++) {
269: PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s.%D\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,j,boffset);
270: boffset += gpiece[r].nvertices*sizeof(PetscScalar) + sizeof(int);
271: }
272: }
273: }
274: PetscFPrintf(PETSC_COMM_SELF,fp," </PointData>\n");
276: PetscFPrintf(PETSC_COMM_SELF,fp," </Piece>\n");
277: }
278: }
280: PetscFPrintf(comm,fp," </UnstructuredGrid>\n");
281: PetscFPrintf(comm,fp," <AppendedData encoding=\"raw\">\n");
282: PetscFPrintf(comm,fp,"_");
284: if (!rank) {
285: PetscInt maxsize = 0;
286: for (r=0; r<size; r++) {
287: maxsize = PetscMax(maxsize, (PetscInt) (gpiece[r].nvertices*3*sizeof(PetscScalar)));
288: maxsize = PetscMax(maxsize, (PetscInt) (gpiece[r].ncells*sizeof(PetscScalar)));
289: maxsize = PetscMax(maxsize, (PetscInt) (gpiece[r].nconn*sizeof(PetscVTKInt)));
290: }
291: PetscMalloc(maxsize,&buffer);
292: }
293: for (r=0; r<size; r++) {
294: if (r == rank) {
295: PetscInt nsend;
296: { /* Position */
297: const PetscScalar *x;
298: PetscScalar *y = NULL;
299: Vec coords;
300: nsend = piece.nvertices*3;
301: DMGetCoordinatesLocal(dm,&coords);
302: VecGetArrayRead(coords,&x);
303: if (dim != 3) {
304: PetscMalloc1(piece.nvertices*3,&y);
305: for (i=0; i<piece.nvertices; i++) {
306: y[i*3+0] = x[i*dim+0];
307: y[i*3+1] = (dim > 1) ? x[i*dim+1] : 0;
308: y[i*3+2] = 0;
309: }
310: }
311: TransferWrite(viewer,fp,r,0,y ? y : x,buffer,nsend,PETSC_SCALAR,tag);
312: PetscFree(y);
313: VecRestoreArrayRead(coords,&x);
314: }
315: { /* Connectivity, offsets, types */
316: PetscVTKInt *connectivity = NULL, *offsets = NULL;
317: PetscVTKType *types = NULL;
318: DMPlexGetVTKConnectivity(dm,&piece,&connectivity,&offsets,&types);
319: TransferWrite(viewer,fp,r,0,connectivity,buffer,piece.nconn,PETSC_INT32,tag);
320: TransferWrite(viewer,fp,r,0,offsets,buffer,piece.ncells,PETSC_INT32,tag);
321: TransferWrite(viewer,fp,r,0,types,buffer,piece.ncells,PETSC_UINT8,tag);
322: PetscFree3(connectivity,offsets,types);
323: }
324: { /* Owners (cell data) */
325: PetscVTKInt *owners;
326: PetscMalloc1(piece.ncells,&owners);
327: for (i=0; i<piece.ncells; i++) owners[i] = rank;
328: TransferWrite(viewer,fp,r,0,owners,buffer,piece.ncells,PETSC_INT32,tag);
329: PetscFree(owners);
330: }
331: /* Cell data */
332: for (link=vtk->link; link; link=link->next) {
333: Vec X = (Vec)link->vec;
334: const PetscScalar *x;
335: PetscScalar *y;
336: PetscInt bs;
337: if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
338: PetscSectionGetDof(dm->defaultSection,cStart,&bs);
339: VecGetArrayRead(X,&x);
340: PetscMalloc1(piece.ncells,&y);
341: for (i=0; i<bs; i++) {
342: PetscInt cnt;
343: for (c=cStart,cnt=0; c<cEnd; c++) {
344: const PetscScalar *xpoint;
345: if (hasLabel) { /* Ignore some cells */
346: PetscInt value;
347: DMPlexGetLabelValue(dm, "vtk", c, &value);
348: if (value != 1) continue;
349: }
350: DMPlexPointLocalRead(dm,c,x,&xpoint);
351: y[cnt++] = xpoint[i];
352: }
353: if (cnt != piece.ncells) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not match");
354: TransferWrite(viewer,fp,r,0,y,buffer,piece.ncells,PETSC_SCALAR,tag);
355: }
356: PetscFree(y);
357: VecRestoreArrayRead(X,&x);
358: }
360: for (link=vtk->link; link; link=link->next) {
361: Vec X = (Vec)link->vec;
362: const PetscScalar *x;
363: PetscScalar *y;
364: PetscInt bs;
365: if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
366: PetscSectionGetDof(dm->defaultSection,vStart,&bs);
367: VecGetArrayRead(X,&x);
368: PetscMalloc1(piece.nvertices,&y);
369: for (i=0; i<bs; i++) {
370: PetscInt cnt;
371: for (v=vStart,cnt=0; v<vEnd; v++) {
372: const PetscScalar *xpoint;
373: DMPlexPointLocalRead(dm,v,x,&xpoint);
374: y[cnt++] = xpoint[i];
375: }
376: if (cnt != piece.nvertices) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not match");
377: TransferWrite(viewer,fp,r,0,y,buffer,piece.nvertices,PETSC_SCALAR,tag);
378: }
379: PetscFree(y);
380: VecRestoreArrayRead(X,&x);
381: }
382: } else if (!rank) {
383: TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].nvertices*3,PETSC_SCALAR,tag); /* positions */
384: TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].nconn,PETSC_INT32,tag); /* connectivity */
385: TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,PETSC_INT32,tag); /* offsets */
386: TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,PETSC_UINT8,tag); /* types */
387: TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,PETSC_INT32,tag); /* owner rank (cells) */
388: /* all cell data */
389: for (link=vtk->link; link; link=link->next) {
390: PetscInt bs;
391: if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
392: PetscSectionGetDof(dm->defaultSection,cStart,&bs);
393: for (i=0; i<bs; i++) {
394: TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,PETSC_SCALAR,tag);
395: }
396: }
397: /* all point data */
398: for (link=vtk->link; link; link=link->next) {
399: PetscInt bs;
400: if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
401: PetscSectionGetDof(dm->defaultSection,vStart,&bs);
402: for (i=0; i<bs; i++) {
403: TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].nvertices,PETSC_SCALAR,tag);
404: }
405: }
406: }
407: }
408: PetscFree(gpiece);
409: PetscFree(buffer);
410: PetscFPrintf(comm,fp,"\n </AppendedData>\n");
411: PetscFPrintf(comm,fp,"</VTKFile>\n");
412: return(0);
413: }