Actual source code: plexvtu.c
petsc-3.13.6 2020-09-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) || defined(PETSC_USE_REAL___FP16)
11: /* output in float if single or half precision in memory */
12: static const char precision[] = "Float32";
13: typedef float PetscVTUReal;
14: #define MPIU_VTUREAL MPI_FLOAT
15: #elif defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
16: /* output in double if double or quad precision in memory */
17: static const char precision[] = "Float64";
18: typedef double PetscVTUReal;
19: #define MPIU_VTUREAL MPI_DOUBLE
20: #else
21: static const char precision[] = "UnknownPrecision";
22: typedef PetscReal PetscVTUReal;
23: #define MPIU_VTUREAL MPIU_REAL
24: #endif
26: static PetscErrorCode TransferWrite(PetscViewer viewer,FILE *fp,PetscMPIInt srank,PetscMPIInt root,const void *send,void *recv,PetscMPIInt count,MPI_Datatype mpidatatype,PetscMPIInt tag)
27: {
28: PetscMPIInt rank;
30: MPI_Comm comm;
33: PetscObjectGetComm((PetscObject)viewer,&comm);
34: MPI_Comm_rank(comm,&rank);
36: if (rank == srank && rank != root) {
37: MPI_Send((void*)send,count,mpidatatype,root,tag,comm);
38: } else if (rank == root) {
39: const void *buffer;
40: if (root == srank) { /* self */
41: buffer = send;
42: } else {
43: MPI_Status status;
44: PetscMPIInt nrecv;
45: MPI_Recv(recv,count,mpidatatype,srank,tag,comm,&status);
46: MPI_Get_count(&status,mpidatatype,&nrecv);
47: if (count != nrecv) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Array size mismatch");
48: buffer = recv;
49: }
50: PetscViewerVTKFWrite(viewer,fp,buffer,count,mpidatatype);
51: }
52: return(0);
53: }
55: static PetscErrorCode DMPlexGetVTKConnectivity(DM dm, PetscBool localized, PieceInfo *piece,PetscVTKInt **oconn,PetscVTKInt **ooffsets,PetscVTKType **otypes)
56: {
58: PetscSection coordSection;
59: PetscVTKInt *conn,*offsets;
60: PetscVTKType *types;
61: PetscInt dim,vStart,vEnd,cStart,cEnd,pStart,pEnd,cellHeight,numLabelCells,hasLabel,c,v,countcell,countconn;
64: PetscMalloc3(piece->nconn,&conn,piece->ncells,&offsets,piece->ncells,&types);
65: DMGetCoordinateSection(dm, &coordSection);
66: DMGetDimension(dm,&dim);
67: DMPlexGetChart(dm,&pStart,&pEnd);
68: DMPlexGetVTKCellHeight(dm, &cellHeight);
69: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
70: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
71: DMGetStratumSize(dm, "vtk", 1, &numLabelCells);
72: hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;
74: countcell = 0;
75: countconn = 0;
76: for (c = cStart; c < cEnd; ++c) {
77: PetscInt nverts,dof = 0,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: if (localized) {
87: PetscSectionGetDof(coordSection, c, &dof);
88: }
89: if (!dof) {
90: PetscInt *closure = NULL;
91: PetscInt closureSize;
93: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
94: for (v = 0; v < closureSize*2; v += 2) {
95: if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
96: if (!localized) conn[countconn++] = closure[v] - vStart;
97: else conn[countconn++] = startoffset + nC;
98: ++nC;
99: }
100: }
101: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
102: } else {
103: for (nC = 0; nC < dof/dim; nC++) conn[countconn++] = startoffset + nC;
104: }
106: {
107: PetscInt n = PetscMin(nC, 8), s = countconn - nC, i, cone[8];
108: for (i = 0; i < n; ++i) cone[i] = conn[s+i];
109: DMPlexReorderCell(dm, c, cone);
110: for (i = 0; i < n; ++i) conn[s+i] = (int)cone[i];
111: }
113: offsets[countcell] = countconn;
115: nverts = countconn - startoffset;
116: DMPlexVTKGetCellType_Internal(dm,dim,nverts,&celltype);
118: types[countcell] = celltype;
119: countcell++;
120: }
121: if (countcell != piece->ncells) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent cell count");
122: if (countconn != piece->nconn) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent connectivity count");
123: *oconn = conn;
124: *ooffsets = offsets;
125: *otypes = types;
126: return(0);
127: }
129: /*
130: Write all fields that have been provided to the viewer
131: Multi-block XML format with binary appended data.
132: */
133: PetscErrorCode DMPlexVTKWriteAll_VTU(DM dm,PetscViewer viewer)
134: {
135: MPI_Comm comm;
136: PetscSection coordSection;
137: PetscViewer_VTK *vtk = (PetscViewer_VTK*)viewer->data;
138: PetscViewerVTKObjectLink link;
139: FILE *fp;
140: PetscMPIInt rank,size,tag;
141: PetscErrorCode ierr;
142: PetscInt dimEmbed,cellHeight,cStart,cEnd,vStart,vEnd,numLabelCells,hasLabel,c,v,r,i;
143: PetscBool localized;
144: PieceInfo piece,*gpiece = NULL;
145: void *buffer = NULL;
146: const char *byte_order = PetscBinaryBigEndian() ? "BigEndian" : "LittleEndian";
147: PetscInt loops_per_scalar;
150: PetscObjectGetComm((PetscObject)dm,&comm);
151: #if defined(PETSC_USE_COMPLEX)
152: loops_per_scalar = 2;
153: #else
154: loops_per_scalar = 1;
155: #endif
156: MPI_Comm_size(comm,&size);
157: MPI_Comm_rank(comm,&rank);
158: tag = ((PetscObject)viewer)->tag;
160: PetscFOpen(comm,vtk->filename,"wb",&fp);
161: PetscFPrintf(comm,fp,"<?xml version=\"1.0\"?>\n");
162: PetscFPrintf(comm,fp,"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"%s\">\n", byte_order);
163: PetscFPrintf(comm,fp," <UnstructuredGrid>\n");
165: DMGetCoordinateDim(dm, &dimEmbed);
166: DMPlexGetVTKCellHeight(dm, &cellHeight);
167: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
168: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
169: DMGetStratumSize(dm, "vtk", 1, &numLabelCells);
170: DMGetCoordinatesLocalized(dm, &localized);
171: DMGetCoordinateSection(dm, &coordSection);
173: hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;
174: piece.nvertices = 0;
175: piece.ncells = 0;
176: piece.nconn = 0;
177: if (!localized) piece.nvertices = vEnd - vStart;
178: for (c = cStart; c < cEnd; ++c) {
179: PetscInt dof = 0;
180: if (hasLabel) {
181: PetscInt value;
183: DMGetLabelValue(dm, "vtk", c, &value);
184: if (value != 1) continue;
185: }
186: if (localized) {
187: PetscSectionGetDof(coordSection, c, &dof);
188: }
189: if (!dof) {
190: PetscInt *closure = NULL;
191: PetscInt closureSize;
193: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
194: for (v = 0; v < closureSize*2; v += 2) {
195: if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
196: piece.nconn++;
197: if (localized) piece.nvertices++;
198: }
199: }
200: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
201: } else {
202: piece.nvertices += dof/dimEmbed;
203: piece.nconn += dof/dimEmbed;
204: }
205: piece.ncells++;
206: }
207: if (!rank) {PetscMalloc1(size,&gpiece);}
208: MPI_Gather((PetscInt*)&piece,sizeof(piece)/sizeof(PetscInt),MPIU_INT,(PetscInt*)gpiece,sizeof(piece)/sizeof(PetscInt),MPIU_INT,0,comm);
210: /*
211: * Write file header
212: */
213: if (!rank) {
214: PetscInt boffset = 0;
216: for (r=0; r<size; r++) {
217: PetscFPrintf(PETSC_COMM_SELF,fp," <Piece NumberOfPoints=\"%D\" NumberOfCells=\"%D\">\n",gpiece[r].nvertices,gpiece[r].ncells);
218: /* Coordinate positions */
219: PetscFPrintf(PETSC_COMM_SELF,fp," <Points>\n");
220: PetscFPrintf(PETSC_COMM_SELF,fp," <DataArray type=\"%s\" Name=\"Position\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,boffset);
221: boffset += gpiece[r].nvertices*3*sizeof(PetscVTUReal) + sizeof(int);
222: PetscFPrintf(PETSC_COMM_SELF,fp," </Points>\n");
223: /* Cell connectivity */
224: PetscFPrintf(PETSC_COMM_SELF,fp," <Cells>\n");
225: PetscFPrintf(PETSC_COMM_SELF,fp," <DataArray type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",boffset);
226: boffset += gpiece[r].nconn*sizeof(PetscInt) + sizeof(int);
227: PetscFPrintf(PETSC_COMM_SELF,fp," <DataArray type=\"Int32\" Name=\"offsets\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",boffset);
228: boffset += gpiece[r].ncells*sizeof(PetscInt) + sizeof(int);
229: PetscFPrintf(PETSC_COMM_SELF,fp," <DataArray type=\"UInt8\" Name=\"types\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",boffset);
230: boffset += gpiece[r].ncells*sizeof(unsigned char) + sizeof(int);
231: PetscFPrintf(PETSC_COMM_SELF,fp," </Cells>\n");
233: /*
234: * Cell Data headers
235: */
236: PetscFPrintf(PETSC_COMM_SELF,fp," <CellData>\n");
237: PetscFPrintf(PETSC_COMM_SELF,fp," <DataArray type=\"Int32\" Name=\"Rank\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",boffset);
238: boffset += gpiece[r].ncells*sizeof(int) + sizeof(int);
239: /* all the vectors */
240: for (link=vtk->link; link; link=link->next) {
241: Vec X = (Vec)link->vec;
242: DM dmX = NULL;
243: PetscInt bs,nfields,field;
244: const char *vecname = "";
245: PetscSection section;
246: if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
247: 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. */
248: PetscObjectGetName((PetscObject)X,&vecname);
249: }
250: VecGetDM(X, &dmX);
251: if (!dmX) dmX = dm;
252: PetscObjectQuery(link->vec, "section", (PetscObject*) §ion);
253: if (!section) {DMGetLocalSection(dmX, §ion);}
254: PetscSectionGetDof(section,cStart,&bs);
255: PetscSectionGetNumFields(section,&nfields);
256: field = 0;
257: if (link->field >= 0) {
258: field = link->field;
259: nfields = field + 1;
260: }
261: for (i=0; field<(nfields?nfields:1); field++) {
262: PetscInt fbs,j;
263: PetscFV fv = NULL;
264: PetscObject f;
265: PetscClassId fClass;
266: const char *fieldname = NULL;
267: char buf[256];
268: PetscBool vector;
269: if (nfields) { /* We have user-defined fields/components */
270: PetscSectionGetFieldDof(section,cStart,field,&fbs);
271: PetscSectionGetFieldName(section,field,&fieldname);
272: } else fbs = bs; /* Say we have one field with 'bs' components */
273: DMGetField(dmX,field,NULL,&f);
274: PetscObjectGetClassId(f,&fClass);
275: if (fClass == PETSCFV_CLASSID) {
276: fv = (PetscFV) f;
277: }
278: if (nfields && !fieldname) {
279: PetscSNPrintf(buf,sizeof(buf),"CellField%D",field);
280: fieldname = buf;
281: }
282: vector = PETSC_FALSE;
283: if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) {
284: vector = PETSC_TRUE;
285: if (fbs > 3) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_ARG_SIZ,"Cell vector fields can have at most 3 components, %D given\n", fbs);
286: for (j = 0; j < fbs; j++) {
287: const char *compName = NULL;
288: if (fv) {
289: PetscFVGetComponentName(fv,j,&compName);
290: if (compName) break;
291: }
292: }
293: if (j < fbs) vector = PETSC_FALSE;
294: }
295: if (vector) {
296: #if defined(PETSC_USE_COMPLEX)
297: PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s.Re\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);
298: boffset += gpiece[r].ncells*3*sizeof(PetscVTUReal) + sizeof(int);
299: PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s.Im\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);
300: boffset += gpiece[r].ncells*3*sizeof(PetscVTUReal) + sizeof(int);
301: #else
302: PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);
303: boffset += gpiece[r].ncells*3*sizeof(PetscVTUReal) + sizeof(int);
304: #endif
305: i+=fbs;
306: } else {
307: for (j=0; j<fbs; j++) {
308: const char *compName = NULL;
309: char finalname[256];
310: if (fv) {
311: PetscFVGetComponentName(fv,j,&compName);
312: }
313: if (compName) {
314: PetscSNPrintf(finalname,255,"%s%s.%s",vecname,fieldname,compName);
315: }
316: else if (fbs > 1) {
317: PetscSNPrintf(finalname,255,"%s%s.%D",vecname,fieldname,j);
318: } else {
319: PetscSNPrintf(finalname,255,"%s%s",vecname,fieldname);
320: }
321: #if defined(PETSC_USE_COMPLEX)
322: PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s.Re\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,finalname,boffset);
323: boffset += gpiece[r].ncells*sizeof(PetscVTUReal) + sizeof(int);
324: PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s.Im\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,finalname,boffset);
325: boffset += gpiece[r].ncells*sizeof(PetscVTUReal) + sizeof(int);
326: #else
327: PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,finalname,boffset);
328: boffset += gpiece[r].ncells*sizeof(PetscVTUReal) + sizeof(int);
329: #endif
330: i++;
331: }
332: }
333: }
334: if (i != bs) SETERRQ2(comm,PETSC_ERR_PLIB,"Total number of field components %D != block size %D",i,bs);
335: }
336: PetscFPrintf(PETSC_COMM_SELF,fp," </CellData>\n");
338: /*
339: * Point Data headers
340: */
341: PetscFPrintf(PETSC_COMM_SELF,fp," <PointData>\n");
342: for (link=vtk->link; link; link=link->next) {
343: Vec X = (Vec)link->vec;
344: DM dmX;
345: PetscInt bs,nfields,field;
346: const char *vecname = "";
347: PetscSection section;
348: if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
349: 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. */
350: PetscObjectGetName((PetscObject)X,&vecname);
351: }
352: VecGetDM(X, &dmX);
353: if (!dmX) dmX = dm;
354: PetscObjectQuery(link->vec, "section", (PetscObject*) §ion);
355: if (!section) {DMGetLocalSection(dmX, §ion);}
356: PetscSectionGetDof(section,vStart,&bs);
357: PetscSectionGetNumFields(section,&nfields);
358: field = 0;
359: if (link->field >= 0) {
360: field = link->field;
361: nfields = field + 1;
362: }
363: for (i=0; field<(nfields?nfields:1); field++) {
364: PetscInt fbs,j;
365: const char *fieldname = NULL;
366: char buf[256];
367: if (nfields) { /* We have user-defined fields/components */
368: PetscSectionGetFieldDof(section,vStart,field,&fbs);
369: PetscSectionGetFieldName(section,field,&fieldname);
370: } else fbs = bs; /* Say we have one field with 'bs' components */
371: if (nfields && !fieldname) {
372: PetscSNPrintf(buf,sizeof(buf),"PointField%D",field);
373: fieldname = buf;
374: }
375: if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) {
376: if (fbs > 3) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_ARG_SIZ,"Point vector fields can have at most 3 components, %D given\n", fbs);
377: #if defined(PETSC_USE_COMPLEX)
378: PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s.Re\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);
379: boffset += gpiece[r].nvertices*3*sizeof(PetscVTUReal) + sizeof(int);
380: PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s.Im\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);
381: boffset += gpiece[r].nvertices*3*sizeof(PetscVTUReal) + sizeof(int);
382: #else
383: PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);
384: boffset += gpiece[r].nvertices*3*sizeof(PetscVTUReal) + sizeof(int);
385: #endif
386: } else {
387: for (j=0; j<fbs; j++) {
388: const char *compName = NULL;
389: char finalname[256];
390: PetscSectionGetComponentName(section,field,j,&compName);
391: PetscSNPrintf(finalname,255,"%s%s.%s",vecname,fieldname,compName);
392: #if defined(PETSC_USE_COMPLEX)
393: PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s.Re\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,finalname,boffset);
394: boffset += gpiece[r].nvertices*sizeof(PetscVTUReal) + sizeof(int);
395: PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s.Im\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,finalname,boffset);
396: boffset += gpiece[r].nvertices*sizeof(PetscVTUReal) + sizeof(int);
397: #else
398: PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,finalname,boffset);
399: boffset += gpiece[r].nvertices*sizeof(PetscVTUReal) + sizeof(int);
400: #endif
401: }
402: }
403: }
404: }
405: PetscFPrintf(PETSC_COMM_SELF,fp," </PointData>\n");
406: PetscFPrintf(PETSC_COMM_SELF,fp," </Piece>\n");
407: }
408: }
410: PetscFPrintf(comm,fp," </UnstructuredGrid>\n");
411: PetscFPrintf(comm,fp," <AppendedData encoding=\"raw\">\n");
412: PetscFPrintf(comm,fp,"_");
414: if (!rank) {
415: PetscInt maxsize = 0;
416: for (r=0; r<size; r++) {
417: maxsize = PetscMax(maxsize, (PetscInt) (gpiece[r].nvertices*3*sizeof(PetscVTUReal)));
418: maxsize = PetscMax(maxsize, (PetscInt) (gpiece[r].ncells*3*sizeof(PetscVTUReal)));
419: maxsize = PetscMax(maxsize, (PetscInt) (gpiece[r].nconn*sizeof(PetscVTKInt)));
420: }
421: PetscMalloc(maxsize,&buffer);
422: }
423: for (r=0; r<size; r++) {
424: if (r == rank) {
425: PetscInt nsend;
426: { /* Position */
427: const PetscScalar *x;
428: PetscVTUReal *y = NULL;
429: Vec coords;
430: PetscBool copy;
432: DMGetCoordinatesLocal(dm,&coords);
433: VecGetArrayRead(coords,&x);
434: #if defined(PETSC_USE_COMPLEX)
435: copy = PETSC_TRUE;
436: #else
437: copy = (PetscBool) (dimEmbed != 3 || localized || (sizeof(PetscReal) != sizeof(PetscVTUReal)));
438: #endif
439: if (copy) {
440: PetscMalloc1(piece.nvertices*3,&y);
441: if (localized) {
442: PetscInt cnt;
443: for (c=cStart,cnt=0; c<cEnd; c++) {
444: PetscInt off, dof;
446: PetscSectionGetDof(coordSection, c, &dof);
447: if (!dof) {
448: PetscInt *closure = NULL;
449: PetscInt closureSize;
451: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
452: for (v = 0; v < closureSize*2; v += 2) {
453: if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
454: PetscSectionGetOffset(coordSection, closure[v], &off);
455: if (dimEmbed != 3) {
456: y[cnt*3+0] = (PetscVTUReal) PetscRealPart(x[off+0]);
457: y[cnt*3+1] = (PetscVTUReal) ((dimEmbed > 1) ? PetscRealPart(x[off+1]) : 0.0);
458: y[cnt*3+2] = (PetscVTUReal) 0.0;
459: } else {
460: y[cnt*3+0] = (PetscVTUReal) PetscRealPart(x[off+0]);
461: y[cnt*3+1] = (PetscVTUReal) PetscRealPart(x[off+1]);
462: y[cnt*3+2] = (PetscVTUReal) PetscRealPart(x[off+2]);
463: }
464: cnt++;
465: }
466: }
467: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
468: } else {
469: PetscSectionGetOffset(coordSection, c, &off);
470: if (dimEmbed != 3) {
471: for (i=0; i<dof/dimEmbed; i++) {
472: y[cnt*3+0] = (PetscVTUReal) PetscRealPart(x[off + i*dimEmbed + 0]);
473: y[cnt*3+1] = (PetscVTUReal) ((dimEmbed > 1) ? PetscRealPart(x[off + i*dimEmbed + 1]) : 0.0);
474: y[cnt*3+2] = (PetscVTUReal) 0.0;
475: cnt++;
476: }
477: } else {
478: for (i=0; i<dof; i ++) {
479: y[cnt*3+i] = (PetscVTUReal) PetscRealPart(x[off + i]);
480: }
481: cnt += dof/dimEmbed;
482: }
483: }
484: }
485: if (cnt != piece.nvertices) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not match");
486: } else {
487: for (i=0; i<piece.nvertices; i++) {
488: y[i*3+0] = (PetscVTUReal) PetscRealPart(x[i*dimEmbed+0]);
489: y[i*3+1] = (PetscVTUReal) ((dimEmbed > 1) ? PetscRealPart(x[i*dimEmbed+1]) : 0.);
490: y[i*3+2] = (PetscVTUReal) ((dimEmbed > 2) ? PetscRealPart(x[i*dimEmbed+2]) : 0.);
491: }
492: }
493: }
494: nsend = piece.nvertices*3;
495: TransferWrite(viewer,fp,r,0,copy ? (const void *) y : (const void *) x,buffer,nsend,MPIU_VTUREAL,tag);
496: PetscFree(y);
497: VecRestoreArrayRead(coords,&x);
498: }
499: { /* Connectivity, offsets, types */
500: PetscVTKInt *connectivity = NULL, *offsets = NULL;
501: PetscVTKType *types = NULL;
502: DMPlexGetVTKConnectivity(dm,localized,&piece,&connectivity,&offsets,&types);
503: TransferWrite(viewer,fp,r,0,connectivity,buffer,piece.nconn,MPI_INT,tag);
504: TransferWrite(viewer,fp,r,0,offsets,buffer,piece.ncells,MPI_INT,tag);
505: TransferWrite(viewer,fp,r,0,types,buffer,piece.ncells,MPI_CHAR,tag);
506: PetscFree3(connectivity,offsets,types);
507: }
508: { /* Owners (cell data) */
509: PetscVTKInt *owners;
510: PetscMalloc1(piece.ncells,&owners);
511: for (i=0; i<piece.ncells; i++) owners[i] = rank;
512: TransferWrite(viewer,fp,r,0,owners,buffer,piece.ncells,MPI_INT,tag);
513: PetscFree(owners);
514: }
515: /* Cell data */
516: for (link=vtk->link; link; link=link->next) {
517: Vec X = (Vec)link->vec;
518: DM dmX;
519: const PetscScalar *x;
520: PetscVTUReal *y;
521: PetscInt bs, nfields, field;
522: PetscSection section = NULL;
524: if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
525: VecGetDM(X, &dmX);
526: if (!dmX) dmX = dm;
527: PetscObjectQuery(link->vec, "section", (PetscObject*) §ion);
528: if (!section) {DMGetLocalSection(dmX, §ion);}
529: PetscSectionGetDof(section,cStart,&bs);
530: PetscSectionGetNumFields(section,&nfields);
531: field = 0;
532: if (link->field >= 0) {
533: field = link->field;
534: nfields = field + 1;
535: }
536: VecGetArrayRead(X,&x);
537: PetscMalloc1(piece.ncells*3,&y);
538: for (i=0; field<(nfields?nfields:1); field++) {
539: PetscInt fbs,j;
540: PetscFV fv = NULL;
541: PetscObject f;
542: PetscClassId fClass;
543: PetscBool vector;
544: if (nfields) { /* We have user-defined fields/components */
545: PetscSectionGetFieldDof(section,cStart,field,&fbs);
546: } else fbs = bs; /* Say we have one field with 'bs' components */
547: DMGetField(dmX,field,NULL,&f);
548: PetscObjectGetClassId(f,&fClass);
549: if (fClass == PETSCFV_CLASSID) {
550: fv = (PetscFV) f;
551: }
552: vector = PETSC_FALSE;
553: if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) {
554: vector = PETSC_TRUE;
555: for (j = 0; j < fbs; j++) {
556: const char *compName = NULL;
557: if (fv) {
558: PetscFVGetComponentName(fv,j,&compName);
559: if (compName) break;
560: }
561: }
562: if (j < fbs) vector = PETSC_FALSE;
563: }
564: if (vector) {
565: PetscInt cnt, l;
566: for (l = 0; l < loops_per_scalar; l++) {
567: for (c=cStart,cnt=0; c<cEnd; c++) {
568: const PetscScalar *xpoint;
569: PetscInt off, j;
571: if (hasLabel) { /* Ignore some cells */
572: PetscInt value;
573: DMGetLabelValue(dmX, "vtk", c, &value);
574: if (value != 1) continue;
575: }
576: if (nfields) {
577: PetscSectionGetFieldOffset(section, c, field, &off);
578: } else {
579: PetscSectionGetOffset(section, c, &off);
580: }
581: xpoint = &x[off];
582: for (j = 0; j < fbs; j++) {
583: y[cnt++] = (PetscVTUReal) (l ? PetscImaginaryPart(xpoint[j]) : PetscRealPart(xpoint[j]));
584: }
585: for (; j < 3; j++) y[cnt++] = 0.;
586: }
587: if (cnt != piece.ncells*3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not match");
588: TransferWrite(viewer,fp,r,0,y,buffer,piece.ncells*3,MPIU_VTUREAL,tag);
589: }
590: } else {
591: for (i=0; i<fbs; i++) {
592: PetscInt cnt, l;
593: for (l = 0; l < loops_per_scalar; l++) {
594: for (c=cStart,cnt=0; c<cEnd; c++) {
595: const PetscScalar *xpoint;
596: PetscInt off;
598: if (hasLabel) { /* Ignore some cells */
599: PetscInt value;
600: DMGetLabelValue(dmX, "vtk", c, &value);
601: if (value != 1) continue;
602: }
603: if (nfields) {
604: PetscSectionGetFieldOffset(section, c, field, &off);
605: } else {
606: PetscSectionGetOffset(section, c, &off);
607: }
608: xpoint = &x[off];
609: y[cnt++] = (PetscVTUReal) (l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i]));
610: }
611: if (cnt != piece.ncells) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not match");
612: TransferWrite(viewer,fp,r,0,y,buffer,piece.ncells,MPIU_VTUREAL,tag);
613: }
614: }
615: }
616: }
617: PetscFree(y);
618: VecRestoreArrayRead(X,&x);
619: }
620: /* point data */
621: for (link=vtk->link; link; link=link->next) {
622: Vec X = (Vec)link->vec;
623: DM dmX;
624: const PetscScalar *x;
625: PetscVTUReal *y;
626: PetscInt bs, nfields, field;
627: PetscSection section = NULL;
629: if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
630: VecGetDM(X, &dmX);
631: if (!dmX) dmX = dm;
632: PetscObjectQuery(link->vec, "section", (PetscObject*) §ion);
633: if (!section) {DMGetLocalSection(dmX, §ion);}
634: PetscSectionGetDof(section,vStart,&bs);
635: PetscSectionGetNumFields(section,&nfields);
636: field = 0;
637: if (link->field >= 0) {
638: field = link->field;
639: nfields = field + 1;
640: }
641: VecGetArrayRead(X,&x);
642: PetscMalloc1(piece.nvertices*3,&y);
643: for (i=0; field<(nfields?nfields:1); field++) {
644: PetscInt fbs,j;
645: if (nfields) { /* We have user-defined fields/components */
646: PetscSectionGetFieldDof(section,vStart,field,&fbs);
647: } else fbs = bs; /* Say we have one field with 'bs' components */
648: if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) {
649: PetscInt cnt, l;
650: for (l = 0; l < loops_per_scalar; l++) {
651: if (!localized) {
652: for (v=vStart,cnt=0; v<vEnd; v++) {
653: PetscInt off;
654: const PetscScalar *xpoint;
656: if (nfields) {
657: PetscSectionGetFieldOffset(section,v,field,&off);
658: } else {
659: PetscSectionGetOffset(section,v,&off);
660: }
661: xpoint = &x[off];
662: for (j = 0; j < fbs; j++) {
663: y[cnt++] = (PetscVTUReal) (l ? PetscImaginaryPart(xpoint[j]) : PetscRealPart(xpoint[j]));
664: }
665: for (; j < 3; j++) y[cnt++] = 0.;
666: }
667: } else {
668: for (c=cStart,cnt=0; c<cEnd; c++) {
669: PetscInt *closure = NULL;
670: PetscInt closureSize, off;
672: DMPlexGetTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure);
673: for (v = 0, off = 0; v < closureSize*2; v += 2) {
674: if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
675: PetscInt voff;
676: const PetscScalar *xpoint;
678: if (nfields) {
679: PetscSectionGetFieldOffset(section,closure[v],field,&voff);
680: } else {
681: PetscSectionGetOffset(section,closure[v],&voff);
682: }
683: xpoint = &x[voff];
684: for (j = 0; j < fbs; j++) {
685: y[cnt + off++] = (PetscVTUReal) (l ? PetscImaginaryPart(xpoint[j]) : PetscRealPart(xpoint[j]));
686: }
687: for (; j < 3; j++) y[cnt + off++] = 0.;
688: }
689: }
690: cnt += off;
691: DMPlexRestoreTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure);
692: }
693: }
694: if (cnt != piece.nvertices*3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not match");
695: TransferWrite(viewer,fp,r,0,y,buffer,piece.nvertices*3,MPIU_VTUREAL,tag);
696: }
697: } else {
698: for (i=0; i<fbs; i++) {
699: PetscInt cnt, l;
700: for (l = 0; l < loops_per_scalar; l++) {
701: if (!localized) {
702: for (v=vStart,cnt=0; v<vEnd; v++) {
703: PetscInt off;
704: const PetscScalar *xpoint;
706: if (nfields) {
707: PetscSectionGetFieldOffset(section,v,field,&off);
708: } else {
709: PetscSectionGetOffset(section,v,&off);
710: }
711: xpoint = &x[off];
712: y[cnt++] = (PetscVTUReal) (l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i]));
713: }
714: } else {
715: for (c=cStart,cnt=0; c<cEnd; c++) {
716: PetscInt *closure = NULL;
717: PetscInt closureSize, off;
719: DMPlexGetTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure);
720: for (v = 0, off = 0; v < closureSize*2; v += 2) {
721: if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
722: PetscInt voff;
723: const PetscScalar *xpoint;
725: if (nfields) {
726: PetscSectionGetFieldOffset(section,closure[v],field,&voff);
727: } else {
728: PetscSectionGetOffset(section,closure[v],&voff);
729: }
730: xpoint = &x[voff];
731: y[cnt + off++] = (l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i]));
732: }
733: }
734: cnt += off;
735: DMPlexRestoreTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure);
736: }
737: }
738: if (cnt != piece.nvertices) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not match");
739: TransferWrite(viewer,fp,r,0,y,buffer,piece.nvertices,MPIU_VTUREAL,tag);
740: }
741: }
742: }
743: }
744: PetscFree(y);
745: VecRestoreArrayRead(X,&x);
746: }
747: } else if (!rank) {
748: PetscInt l;
750: TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].nvertices*3,MPIU_VTUREAL,tag); /* positions */
751: TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].nconn,MPI_INT,tag); /* connectivity */
752: TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,MPI_INT,tag); /* offsets */
753: TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,MPI_CHAR,tag); /* types */
754: TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,MPI_INT,tag); /* owner rank (cells) */
755: /* all cell data */
756: for (link=vtk->link; link; link=link->next) {
757: Vec X = (Vec)link->vec;
758: PetscInt bs, nfields, field;
759: DM dmX;
760: PetscSection section = NULL;
762: if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
763: VecGetDM(X, &dmX);
764: if (!dmX) dmX = dm;
765: PetscObjectQuery(link->vec, "section", (PetscObject*) §ion);
766: if (!section) {DMGetLocalSection(dmX, §ion);}
767: PetscSectionGetDof(section,cStart,&bs);
768: PetscSectionGetNumFields(section,&nfields);
769: field = 0;
770: if (link->field >= 0) {
771: field = link->field;
772: nfields = field + 1;
773: }
774: for (i=0; field<(nfields?nfields:1); field++) {
775: PetscInt fbs,j;
776: PetscFV fv = NULL;
777: PetscObject f;
778: PetscClassId fClass;
779: PetscBool vector;
780: if (nfields) { /* We have user-defined fields/components */
781: PetscSectionGetFieldDof(section,cStart,field,&fbs);
782: } else fbs = bs; /* Say we have one field with 'bs' components */
783: DMGetField(dmX,field,NULL,&f);
784: PetscObjectGetClassId(f,&fClass);
785: if (fClass == PETSCFV_CLASSID) {
786: fv = (PetscFV) f;
787: }
788: vector = PETSC_FALSE;
789: if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) {
790: vector = PETSC_TRUE;
791: for (j = 0; j < fbs; j++) {
792: const char *compName = NULL;
793: if (fv) {
794: PetscFVGetComponentName(fv,j,&compName);
795: if (compName) break;
796: }
797: }
798: if (j < fbs) vector = PETSC_FALSE;
799: }
800: if (vector) {
801: for (l = 0; l < loops_per_scalar; l++) {
802: TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells*3,MPIU_VTUREAL,tag);
803: }
804: } else {
805: for (i=0; i<fbs; i++) {
806: for (l = 0; l < loops_per_scalar; l++) {
807: TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,MPIU_VTUREAL,tag);
808: }
809: }
810: }
811: }
812: }
813: /* all point data */
814: for (link=vtk->link; link; link=link->next) {
815: Vec X = (Vec)link->vec;
816: DM dmX;
817: PetscInt bs, nfields, field;
818: PetscSection section = NULL;
820: if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
821: VecGetDM(X, &dmX);
822: if (!dmX) dmX = dm;
823: PetscObjectQuery(link->vec, "section", (PetscObject*) §ion);
824: if (!section) {DMGetLocalSection(dmX, §ion);}
825: PetscSectionGetDof(section,vStart,&bs);
826: PetscSectionGetNumFields(section,&nfields);
827: field = 0;
828: if (link->field >= 0) {
829: field = link->field;
830: nfields = field + 1;
831: }
832: for (i=0; field<(nfields?nfields:1); field++) {
833: PetscInt fbs;
834: if (nfields) { /* We have user-defined fields/components */
835: PetscSectionGetFieldDof(section,vStart,field,&fbs);
836: } else fbs = bs; /* Say we have one field with 'bs' components */
837: if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) {
838: for (l = 0; l < loops_per_scalar; l++) {
839: TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].nvertices*3,MPIU_VTUREAL,tag);
840: }
841: } else {
842: for (i=0; i<fbs; i++) {
843: for (l = 0; l < loops_per_scalar; l++) {
844: TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].nvertices,MPIU_VTUREAL,tag);
845: }
846: }
847: }
848: }
849: }
850: }
851: }
852: PetscFree(gpiece);
853: PetscFree(buffer);
854: PetscFPrintf(comm,fp,"\n </AppendedData>\n");
855: PetscFPrintf(comm,fp,"</VTKFile>\n");
856: PetscFClose(comm,fp);
857: return(0);
858: }