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) PetscCall(PetscSectionGetFieldOffset(section, c, field, &off));
591: else PetscCall(PetscSectionGetOffset(section, c, &off));
592: xpoint = &x[off];
593: for (j = 0; j < fbs; j++) y[cnt++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[j]) : PetscRealPart(xpoint[j]));
594: for (; j < 3; j++) y[cnt++] = 0.;
595: }
596: PetscCheck(cnt == piece.ncells * 3, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count does not match");
597: PetscCall(TransferWrite(comm, viewer, fp, r, 0, y, buffer, piece.ncells * 3, MPIU_VTUREAL, tag));
598: }
599: } else {
600: for (i = 0; i < fbs; i++) {
601: PetscInt cnt, l;
602: for (l = 0; l < loops_per_scalar; l++) {
603: for (c = cStart, cnt = 0; c < cEnd; c++) {
604: const PetscScalar *xpoint;
605: PetscInt off;
607: if (hasLabel) { /* Ignore some cells */
608: PetscInt value;
609: PetscCall(DMGetLabelValue(dmX, "vtk", c, &value));
610: if (value != 1) continue;
611: }
612: if (nfields) PetscCall(PetscSectionGetFieldOffset(section, c, field, &off));
613: else PetscCall(PetscSectionGetOffset(section, c, &off));
614: xpoint = &x[off];
615: y[cnt++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i]));
616: }
617: PetscCheck(cnt == piece.ncells, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count does not match");
618: PetscCall(TransferWrite(comm, viewer, fp, r, 0, y, buffer, piece.ncells, MPIU_VTUREAL, tag));
619: }
620: }
621: }
622: }
623: PetscCall(PetscFree(y));
624: PetscCall(VecRestoreArrayRead(X, &x));
625: }
626: /* point data */
627: for (link = vtk->link; link; link = link->next) {
628: Vec X = (Vec)link->vec;
629: DM dmX;
630: const PetscScalar *x;
631: PetscVTUReal *y;
632: PetscInt bs = 1, nfields, field;
633: PetscSection section = NULL;
635: if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
636: PetscCall(VecGetDM(X, &dmX));
637: if (!dmX) dmX = dm;
638: PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)§ion));
639: if (!section) PetscCall(DMGetLocalSection(dmX, §ion));
640: if (vEnd > vStart) PetscCall(PetscSectionGetDof(section, vStart, &bs));
641: PetscCall(PetscSectionGetNumFields(section, &nfields));
642: field = 0;
643: if (link->field >= 0) {
644: field = link->field;
645: nfields = field + 1;
646: }
647: PetscCall(VecGetArrayRead(X, &x));
648: PetscCall(PetscMalloc1(piece.nvertices * 3, &y));
649: for (; field < (nfields ? nfields : 1); field++) {
650: PetscInt fbs, j;
651: if (nfields && vEnd > vStart) { /* We have user-defined fields/components */
652: PetscCall(PetscSectionGetFieldDof(section, vStart, field, &fbs));
653: } else fbs = bs; /* Say we have one field with 'bs' components */
654: if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) {
655: PetscInt cnt, l;
656: for (l = 0; l < loops_per_scalar; l++) {
657: if (!localized) {
658: for (v = vStart, cnt = 0; v < vEnd; v++) {
659: PetscInt off;
660: const PetscScalar *xpoint;
662: if (nfields) PetscCall(PetscSectionGetFieldOffset(section, v, field, &off));
663: else PetscCall(PetscSectionGetOffset(section, v, &off));
664: xpoint = &x[off];
665: for (j = 0; j < fbs; j++) y[cnt++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[j]) : PetscRealPart(xpoint[j]));
666: for (; j < 3; j++) y[cnt++] = 0.;
667: }
668: } else {
669: for (c = cStart, cnt = 0; c < cEnd; c++) {
670: PetscInt *closure = NULL;
671: PetscInt closureSize, off;
673: PetscCall(DMPlexGetTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure));
674: for (v = 0, off = 0; v < closureSize * 2; v += 2) {
675: if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
676: PetscInt voff;
677: const PetscScalar *xpoint;
679: if (nfields) {
680: PetscCall(PetscSectionGetFieldOffset(section, closure[v], field, &voff));
681: } else {
682: PetscCall(PetscSectionGetOffset(section, closure[v], &voff));
683: }
684: xpoint = &x[voff];
685: for (j = 0; j < fbs; j++) y[cnt + off++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[j]) : PetscRealPart(xpoint[j]));
686: for (; j < 3; j++) y[cnt + off++] = 0.;
687: }
688: }
689: cnt += off;
690: PetscCall(DMPlexRestoreTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure));
691: }
692: }
693: PetscCheck(cnt == piece.nvertices * 3, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count does not match");
694: PetscCall(TransferWrite(comm, viewer, fp, r, 0, y, buffer, piece.nvertices * 3, MPIU_VTUREAL, tag));
695: }
696: } else {
697: for (i = 0; i < fbs; i++) {
698: PetscInt cnt, l;
699: for (l = 0; l < loops_per_scalar; l++) {
700: if (!localized) {
701: for (v = vStart, cnt = 0; v < vEnd; v++) {
702: PetscInt off;
703: const PetscScalar *xpoint;
705: if (nfields) PetscCall(PetscSectionGetFieldOffset(section, v, field, &off));
706: else PetscCall(PetscSectionGetOffset(section, v, &off));
707: xpoint = &x[off];
708: y[cnt++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i]));
709: }
710: } else {
711: for (c = cStart, cnt = 0; c < cEnd; c++) {
712: PetscInt *closure = NULL;
713: PetscInt closureSize, off;
715: PetscCall(DMPlexGetTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure));
716: for (v = 0, off = 0; v < closureSize * 2; v += 2) {
717: if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
718: PetscInt voff;
719: const PetscScalar *xpoint;
721: if (nfields) PetscCall(PetscSectionGetFieldOffset(section, closure[v], field, &voff));
722: else PetscCall(PetscSectionGetOffset(section, closure[v], &voff));
723: xpoint = &x[voff];
724: y[cnt + off++] = (l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i]));
725: }
726: }
727: cnt += off;
728: PetscCall(DMPlexRestoreTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure));
729: }
730: }
731: PetscCheck(cnt == piece.nvertices, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count does not match");
732: PetscCall(TransferWrite(comm, viewer, fp, r, 0, y, buffer, piece.nvertices, MPIU_VTUREAL, tag));
733: }
734: }
735: }
736: }
737: PetscCall(PetscFree(y));
738: PetscCall(VecRestoreArrayRead(X, &x));
739: }
740: } else if (rank == 0) {
741: PetscInt l;
743: PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].nvertices * 3, MPIU_VTUREAL, tag)); /* positions */
744: PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].nconn, MPI_INT, tag)); /* connectivity */
745: PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells, MPI_INT, tag)); /* offsets */
746: PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells, MPI_CHAR, tag)); /* types */
747: PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells, MPI_INT, tag)); /* owner rank (cells) */
748: /* all cell data */
749: for (link = vtk->link; link; link = link->next) {
750: Vec X = (Vec)link->vec;
751: PetscInt bs = 1, nfields, field;
752: DM dmX;
753: PetscSection section = NULL;
755: if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
756: PetscCall(VecGetDM(X, &dmX));
757: if (!dmX) dmX = dm;
758: PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)§ion));
759: if (!section) PetscCall(DMGetLocalSection(dmX, §ion));
760: if (cEnd > cStart) PetscCall(PetscSectionGetDof(section, cStart, &bs));
761: PetscCall(PetscSectionGetNumFields(section, &nfields));
762: field = 0;
763: if (link->field >= 0) {
764: field = link->field;
765: nfields = field + 1;
766: }
767: for (; field < (nfields ? nfields : 1); field++) {
768: PetscInt fbs, j;
769: PetscFV fv = NULL;
770: PetscBool vector;
772: if (nfields && cEnd > cStart) { /* We have user-defined fields/components */
773: PetscCall(PetscSectionGetFieldDof(section, cStart, field, &fbs));
774: } else fbs = bs; /* Say we have one field with 'bs' components */
775: PetscCall(DMGetFieldIfFV_Private(dmX, field, &fv));
776: vector = PETSC_FALSE;
777: if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) {
778: vector = PETSC_TRUE;
779: for (j = 0; j < fbs; j++) {
780: const char *compName = NULL;
781: if (fv) {
782: PetscCall(PetscFVGetComponentName(fv, j, &compName));
783: if (compName) break;
784: }
785: }
786: if (j < fbs) vector = PETSC_FALSE;
787: }
788: if (vector) {
789: for (l = 0; l < loops_per_scalar; l++) PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells * 3, MPIU_VTUREAL, tag));
790: } else {
791: for (i = 0; i < fbs; i++) {
792: for (l = 0; l < loops_per_scalar; l++) PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells, MPIU_VTUREAL, tag));
793: }
794: }
795: }
796: }
797: /* all point data */
798: for (link = vtk->link; link; link = link->next) {
799: Vec X = (Vec)link->vec;
800: DM dmX;
801: PetscInt bs = 1, nfields, field;
802: PetscSection section = NULL;
804: if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
805: PetscCall(VecGetDM(X, &dmX));
806: if (!dmX) dmX = dm;
807: PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)§ion));
808: if (!section) PetscCall(DMGetLocalSection(dmX, §ion));
809: if (vEnd > vStart) PetscCall(PetscSectionGetDof(section, vStart, &bs));
810: PetscCall(PetscSectionGetNumFields(section, &nfields));
811: field = 0;
812: if (link->field >= 0) {
813: field = link->field;
814: nfields = field + 1;
815: }
816: for (; field < (nfields ? nfields : 1); field++) {
817: PetscInt fbs;
818: if (nfields && vEnd > vStart) { /* We have user-defined fields/components */
819: PetscCall(PetscSectionGetFieldDof(section, vStart, field, &fbs));
820: } else fbs = bs; /* Say we have one field with 'bs' components */
821: if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) {
822: for (l = 0; l < loops_per_scalar; l++) PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].nvertices * 3, MPIU_VTUREAL, tag));
823: } else {
824: for (i = 0; i < fbs; i++) {
825: for (l = 0; l < loops_per_scalar; l++) PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].nvertices, MPIU_VTUREAL, tag));
826: }
827: }
828: }
829: }
830: }
831: }
832: PetscCall(PetscFree(gpiece));
833: PetscCall(PetscFree(buffer));
834: PetscCall(PetscFPrintf(comm, fp, "\n </AppendedData>\n"));
835: PetscCall(PetscFPrintf(comm, fp, "</VTKFile>\n"));
836: PetscCall(PetscFClose(comm, fp));
837: finalize:
838: /* this code sends to rank 0 that writes.
839: It may lead to very unbalanced log_view timings
840: of the next PETSc function logged.
841: Since this call is not performance critical, we
842: issue a barrier here to synchronize the processes */
843: PetscCall(PetscBarrier((PetscObject)viewer));
844: PetscFunctionReturn(PETSC_SUCCESS);
845: }