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