Actual source code: plexvtu.c
petsc-3.7.3 2016-08-01
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: DMGetDimension(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: DMGetStratumSize(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: DMGetLabelValue(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 dimEmbed,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: DMGetCoordinateDim(dm, &dimEmbed);
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: DMGetStratumSize(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: DMGetLabelValue(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: PetscFV fv = NULL;
223: PetscObject f;
224: PetscClassId fClass;
225: const char *fieldname = NULL;
226: char buf[256];
227: if (nfields) { /* We have user-defined fields/components */
228: PetscSectionGetFieldDof(dm->defaultSection,cStart,field,&fbs);
229: PetscSectionGetFieldName(dm->defaultSection,field,&fieldname);
230: } else fbs = bs; /* Say we have one field with 'bs' components */
231: DMGetField(dm,field,&f);
232: PetscObjectGetClassId(f,&fClass);
233: if (fClass == PETSCFV_CLASSID) {
234: fv = (PetscFV) f;
235: }
236: if (!fieldname) {
237: PetscSNPrintf(buf,sizeof(buf),"CellField%D",field);
238: fieldname = buf;
239: }
240: for (j=0; j<fbs; j++) {
241: const char *compName = NULL;
242: if (fv) {
243: PetscFVGetComponentName(fv,j,&compName);
244: }
245: if (compName) {
246: PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s.%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,compName,boffset);
247: }
248: else {
249: PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s.%D\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,j,boffset);
250: }
251: boffset += gpiece[r].ncells*sizeof(PetscScalar) + sizeof(int);
252: i++;
253: }
254: }
255: if (i != bs) SETERRQ2(comm,PETSC_ERR_PLIB,"Total number of field components %D != block size %D",i,bs);
256: }
257: PetscFPrintf(PETSC_COMM_SELF,fp," </CellData>\n");
259: /*
260: * Point Data headers
261: */
262: PetscFPrintf(PETSC_COMM_SELF,fp," <PointData>\n");
263: for (link=vtk->link; link; link=link->next) {
264: Vec X = (Vec)link->vec;
265: PetscInt bs,nfields,field;
266: const char *vecname = "";
267: if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
268: 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. */
269: PetscObjectGetName((PetscObject)X,&vecname);
270: }
271: PetscSectionGetDof(dm->defaultSection,vStart,&bs);
272: PetscSectionGetNumFields(dm->defaultSection,&nfields);
273: for (field=0,i=0; field<(nfields?nfields:1); field++) {
274: PetscInt fbs,j;
275: const char *fieldname = NULL;
276: char buf[256];
277: if (nfields) { /* We have user-defined fields/components */
278: PetscSectionGetFieldDof(dm->defaultSection,vStart,field,&fbs);
279: PetscSectionGetFieldName(dm->defaultSection,field,&fieldname);
280: } else fbs = bs; /* Say we have one field with 'bs' components */
281: if (!fieldname) {
282: PetscSNPrintf(buf,sizeof(buf),"PointField%D",field);
283: fieldname = buf;
284: }
285: for (j=0; j<fbs; j++) {
286: PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s.%D\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,j,boffset);
287: boffset += gpiece[r].nvertices*sizeof(PetscScalar) + sizeof(int);
288: }
289: }
290: }
291: PetscFPrintf(PETSC_COMM_SELF,fp," </PointData>\n");
293: PetscFPrintf(PETSC_COMM_SELF,fp," </Piece>\n");
294: }
295: }
297: PetscFPrintf(comm,fp," </UnstructuredGrid>\n");
298: PetscFPrintf(comm,fp," <AppendedData encoding=\"raw\">\n");
299: PetscFPrintf(comm,fp,"_");
301: if (!rank) {
302: PetscInt maxsize = 0;
303: for (r=0; r<size; r++) {
304: maxsize = PetscMax(maxsize, (PetscInt) (gpiece[r].nvertices*3*sizeof(PetscScalar)));
305: maxsize = PetscMax(maxsize, (PetscInt) (gpiece[r].ncells*sizeof(PetscScalar)));
306: maxsize = PetscMax(maxsize, (PetscInt) (gpiece[r].nconn*sizeof(PetscVTKInt)));
307: }
308: PetscMalloc(maxsize,&buffer);
309: }
310: for (r=0; r<size; r++) {
311: if (r == rank) {
312: PetscInt nsend;
313: { /* Position */
314: const PetscScalar *x;
315: PetscScalar *y = NULL;
316: Vec coords;
317: nsend = piece.nvertices*3;
318: DMGetCoordinatesLocal(dm,&coords);
319: VecGetArrayRead(coords,&x);
320: if (dimEmbed != 3) {
321: PetscMalloc1(piece.nvertices*3,&y);
322: for (i=0; i<piece.nvertices; i++) {
323: y[i*3+0] = x[i*dimEmbed+0];
324: y[i*3+1] = (dimEmbed > 1) ? x[i*dimEmbed+1] : 0;
325: y[i*3+2] = 0;
326: }
327: }
328: TransferWrite(viewer,fp,r,0,y ? y : x,buffer,nsend,PETSC_SCALAR,tag);
329: PetscFree(y);
330: VecRestoreArrayRead(coords,&x);
331: }
332: { /* Connectivity, offsets, types */
333: PetscVTKInt *connectivity = NULL, *offsets = NULL;
334: PetscVTKType *types = NULL;
335: DMPlexGetVTKConnectivity(dm,&piece,&connectivity,&offsets,&types);
336: TransferWrite(viewer,fp,r,0,connectivity,buffer,piece.nconn,PETSC_INT32,tag);
337: TransferWrite(viewer,fp,r,0,offsets,buffer,piece.ncells,PETSC_INT32,tag);
338: TransferWrite(viewer,fp,r,0,types,buffer,piece.ncells,PETSC_UINT8,tag);
339: PetscFree3(connectivity,offsets,types);
340: }
341: { /* Owners (cell data) */
342: PetscVTKInt *owners;
343: PetscMalloc1(piece.ncells,&owners);
344: for (i=0; i<piece.ncells; i++) owners[i] = rank;
345: TransferWrite(viewer,fp,r,0,owners,buffer,piece.ncells,PETSC_INT32,tag);
346: PetscFree(owners);
347: }
348: /* Cell data */
349: for (link=vtk->link; link; link=link->next) {
350: Vec X = (Vec)link->vec;
351: const PetscScalar *x;
352: PetscScalar *y;
353: PetscInt bs;
354: if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
355: PetscSectionGetDof(dm->defaultSection,cStart,&bs);
356: VecGetArrayRead(X,&x);
357: PetscMalloc1(piece.ncells,&y);
358: for (i=0; i<bs; i++) {
359: PetscInt cnt;
360: for (c=cStart,cnt=0; c<cEnd; c++) {
361: PetscScalar *xpoint;
362: if (hasLabel) { /* Ignore some cells */
363: PetscInt value;
364: DMGetLabelValue(dm, "vtk", c, &value);
365: if (value != 1) continue;
366: }
367: DMPlexPointLocalRead(dm,c,x,&xpoint);
368: y[cnt++] = xpoint[i];
369: }
370: if (cnt != piece.ncells) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not match");
371: TransferWrite(viewer,fp,r,0,y,buffer,piece.ncells,PETSC_SCALAR,tag);
372: }
373: PetscFree(y);
374: VecRestoreArrayRead(X,&x);
375: }
377: for (link=vtk->link; link; link=link->next) {
378: Vec X = (Vec)link->vec;
379: const PetscScalar *x;
380: PetscScalar *y;
381: PetscInt bs;
382: if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
383: PetscSectionGetDof(dm->defaultSection,vStart,&bs);
384: VecGetArrayRead(X,&x);
385: PetscMalloc1(piece.nvertices,&y);
386: for (i=0; i<bs; i++) {
387: PetscInt cnt;
388: for (v=vStart,cnt=0; v<vEnd; v++) {
389: PetscScalar *xpoint;
390: DMPlexPointLocalRead(dm,v,x,&xpoint);
391: y[cnt++] = xpoint[i];
392: }
393: if (cnt != piece.nvertices) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not match");
394: TransferWrite(viewer,fp,r,0,y,buffer,piece.nvertices,PETSC_SCALAR,tag);
395: }
396: PetscFree(y);
397: VecRestoreArrayRead(X,&x);
398: }
399: } else if (!rank) {
400: TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].nvertices*3,PETSC_SCALAR,tag); /* positions */
401: TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].nconn,PETSC_INT32,tag); /* connectivity */
402: TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,PETSC_INT32,tag); /* offsets */
403: TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,PETSC_UINT8,tag); /* types */
404: TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,PETSC_INT32,tag); /* owner rank (cells) */
405: /* all cell data */
406: for (link=vtk->link; link; link=link->next) {
407: PetscInt bs;
408: if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
409: PetscSectionGetDof(dm->defaultSection,cStart,&bs);
410: for (i=0; i<bs; i++) {
411: TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,PETSC_SCALAR,tag);
412: }
413: }
414: /* all point data */
415: for (link=vtk->link; link; link=link->next) {
416: PetscInt bs;
417: if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
418: PetscSectionGetDof(dm->defaultSection,vStart,&bs);
419: for (i=0; i<bs; i++) {
420: TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].nvertices,PETSC_SCALAR,tag);
421: }
422: }
423: }
424: }
425: PetscFree(gpiece);
426: PetscFree(buffer);
427: PetscFPrintf(comm,fp,"\n </AppendedData>\n");
428: PetscFPrintf(comm,fp,"</VTKFile>\n");
429: PetscFClose(comm,fp);
430: return(0);
431: }