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, PetscMPIInt 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(MPI_Send((void *)send, count, mpidatatype, root, tag, comm));
 34:   } else if (rank == root) {
 35:     const void *buffer;
 36:     if (root == srank) { /* self */
 37:       buffer = send;
 38:     } else {
 39:       MPI_Status  status;
 40:       PetscMPIInt nrecv;
 41:       PetscCallMPI(MPI_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) conn[countconn++] = closure[v] - vStart;
 91:           else conn[countconn++] = startoffset + nC;
 92:           ++nC;
 93:         }
 94:       }
 95:       PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
 96:     } else {
 97:       for (nC = 0; nC < dof / dim; nC++) conn[countconn++] = startoffset + nC;
 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) conn[s + i] = (int)cone[i];
105:     }

107:     offsets[countcell] = countconn;

109:     nverts = countconn - startoffset;
110:     PetscCall(DMPlexVTKGetCellType_Internal(dm, dim, nverts, &celltype));

112:     types[countcell] = celltype;
113:     countcell++;
114:   }
115:   PetscCheck(countcell == piece->ncells, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent cell count");
116:   PetscCheck(countconn == piece->nconn, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent connectivity count");
117:   *oconn    = conn;
118:   *ooffsets = offsets;
119:   *otypes   = types;
120:   PetscFunctionReturn(PETSC_SUCCESS);
121: }

123: /*
124:   Write all fields that have been provided to the viewer
125:   Multi-block XML format with binary appended data.
126: */
127: PetscErrorCode DMPlexVTKWriteAll_VTU(DM dm, PetscViewer viewer)
128: {
129:   MPI_Comm                 comm;
130:   PetscSection             coordSection, cellCoordSection;
131:   PetscViewer_VTK         *vtk = (PetscViewer_VTK *)viewer->data;
132:   PetscViewerVTKObjectLink link;
133:   FILE                    *fp;
134:   PetscMPIInt              rank, size, tag;
135:   PetscInt                 dimEmbed, cellHeight, cStart, cEnd, vStart, vEnd, numLabelCells, hasLabel, c, v, r, i;
136:   PetscBool                localized;
137:   PieceInfo                piece, *gpiece = NULL;
138:   void                    *buffer     = NULL;
139:   const char              *byte_order = PetscBinaryBigEndian() ? "BigEndian" : "LittleEndian";
140:   PetscInt                 loops_per_scalar;

142:   PetscFunctionBegin;
143:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
144: #if defined(PETSC_USE_COMPLEX)
145:   loops_per_scalar = 2;
146: #else
147:   loops_per_scalar = 1;
148: #endif
149:   PetscCallMPI(MPI_Comm_size(comm, &size));
150:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
151:   PetscCall(PetscCommGetNewTag(comm, &tag));

153:   PetscCall(PetscFOpen(comm, vtk->filename, "wb", &fp));
154:   PetscCall(PetscFPrintf(comm, fp, "<?xml version=\"1.0\"?>\n"));
155:   PetscCall(PetscFPrintf(comm, fp, "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"%s\" header_type=\"UInt64\">\n", byte_order));
156:   PetscCall(PetscFPrintf(comm, fp, "  <UnstructuredGrid>\n"));

158:   PetscCall(DMGetCoordinateDim(dm, &dimEmbed));
159:   PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
160:   PetscCall(DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd));
161:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
162:   PetscCall(DMGetStratumSize(dm, "vtk", 1, &numLabelCells));
163:   PetscCall(DMGetCoordinatesLocalized(dm, &localized));
164:   PetscCall(DMGetCoordinateSection(dm, &coordSection));
165:   PetscCall(DMGetCellCoordinateSection(dm, &cellCoordSection));

167:   hasLabel        = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;
168:   piece.nvertices = 0;
169:   piece.ncells    = 0;
170:   piece.nconn     = 0;
171:   if (!localized) piece.nvertices = vEnd - vStart;
172:   for (c = cStart; c < cEnd; ++c) {
173:     PetscInt dof = 0;
174:     if (hasLabel) {
175:       PetscInt value;

177:       PetscCall(DMGetLabelValue(dm, "vtk", c, &value));
178:       if (value != 1) continue;
179:     }
180:     if (localized) PetscCall(PetscSectionGetDof(cellCoordSection, c, &dof));
181:     if (!dof) {
182:       PetscInt *closure = NULL;
183:       PetscInt  closureSize;

185:       PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
186:       for (v = 0; v < closureSize * 2; v += 2) {
187:         if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
188:           piece.nconn++;
189:           if (localized) piece.nvertices++;
190:         }
191:       }
192:       PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
193:     } else {
194:       piece.nvertices += dof / dimEmbed;
195:       piece.nconn += dof / dimEmbed;
196:     }
197:     piece.ncells++;
198:   }
199:   if (rank == 0) PetscCall(PetscMalloc1(size, &gpiece));
200:   PetscCallMPI(MPI_Gather((PetscInt *)&piece, sizeof(piece) / sizeof(PetscInt), MPIU_INT, (PetscInt *)gpiece, sizeof(piece) / sizeof(PetscInt), MPIU_INT, 0, comm));

202:   /*
203:    * Write file header
204:    */
205:   if (rank == 0) {
206:     PetscInt64 boffset = 0;

208:     for (r = 0; r < size; r++) {
209:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "    <Piece NumberOfPoints=\"%" PetscInt_FMT "\" NumberOfCells=\"%" PetscInt_FMT "\">\n", gpiece[r].nvertices, gpiece[r].ncells));
210:       /* Coordinate positions */
211:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "      <Points>\n"));
212:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "        <DataArray type=\"%s\" Name=\"Position\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, boffset));
213:       boffset += gpiece[r].nvertices * 3 * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0);
214:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "      </Points>\n"));
215:       /* Cell connectivity */
216:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "      <Cells>\n"));
217:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "        <DataArray type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", boffset));
218:       boffset += gpiece[r].nconn * sizeof(PetscVTKInt) + (gpiece[r].nconn ? sizeof(PetscInt64) : 0);
219:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "        <DataArray type=\"Int32\" Name=\"offsets\"      NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", boffset));
220:       boffset += gpiece[r].ncells * sizeof(PetscVTKInt) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0);
221:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "        <DataArray type=\"UInt8\" Name=\"types\"        NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", boffset));
222:       boffset += gpiece[r].ncells * sizeof(unsigned char) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0);
223:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "      </Cells>\n"));

225:       /*
226:        * Cell Data headers
227:        */
228:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "      <CellData>\n"));
229:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "        <DataArray type=\"Int32\" Name=\"Rank\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", boffset));
230:       boffset += gpiece[r].ncells * sizeof(PetscVTKInt) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0);
231:       /* all the vectors */
232:       for (link = vtk->link; link; link = link->next) {
233:         Vec          X       = (Vec)link->vec;
234:         DM           dmX     = NULL;
235:         PetscInt     bs      = 1, nfields, field;
236:         const char  *vecname = "";
237:         PetscSection section;
238:         if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
239:         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. */
240:           PetscCall(PetscObjectGetName((PetscObject)X, &vecname));
241:         }
242:         PetscCall(VecGetDM(X, &dmX));
243:         if (!dmX) dmX = dm;
244:         PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)&section));
245:         if (!section) PetscCall(DMGetLocalSection(dmX, &section));
246:         if (cEnd > cStart) PetscCall(PetscSectionGetDof(section, cStart, &bs));
247:         PetscCall(PetscSectionGetNumFields(section, &nfields));
248:         field = 0;
249:         if (link->field >= 0) {
250:           field   = link->field;
251:           nfields = field + 1;
252:         }
253:         for (i = 0; field < (nfields ? nfields : 1); field++) {
254:           PetscInt     fbs, j;
255:           PetscFV      fv = NULL;
256:           PetscObject  f;
257:           PetscClassId fClass;
258:           const char  *fieldname = NULL;
259:           char         buf[256];
260:           PetscBool    vector;
261:           if (nfields) { /* We have user-defined fields/components */
262:             PetscCall(PetscSectionGetFieldDof(section, cStart, field, &fbs));
263:             PetscCall(PetscSectionGetFieldName(section, field, &fieldname));
264:           } else fbs = bs; /* Say we have one field with 'bs' components */
265:           PetscCall(DMGetField(dmX, field, NULL, &f));
266:           PetscCall(PetscObjectGetClassId(f, &fClass));
267:           if (fClass == PETSCFV_CLASSID) fv = (PetscFV)f;
268:           if (nfields && !fieldname) {
269:             PetscCall(PetscSNPrintf(buf, sizeof(buf), "CellField%" PetscInt_FMT, field));
270:             fieldname = buf;
271:           }
272:           vector = PETSC_FALSE;
273:           if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) {
274:             vector = PETSC_TRUE;
275:             PetscCheck(fbs <= 3, PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_SIZ, "Cell vector fields can have at most 3 components, %" PetscInt_FMT " given", fbs);
276:             for (j = 0; j < fbs; j++) {
277:               const char *compName = NULL;
278:               if (fv) {
279:                 PetscCall(PetscFVGetComponentName(fv, j, &compName));
280:                 if (compName) break;
281:               }
282:             }
283:             if (j < fbs) vector = PETSC_FALSE;
284:           }
285:           if (vector) {
286: #if defined(PETSC_USE_COMPLEX)
287:             PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s%s.Re\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset));
288:             boffset += gpiece[r].ncells * 3 * sizeof(PetscVTUReal) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0);
289:             PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s%s.Im\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset));
290:             boffset += gpiece[r].ncells * 3 * sizeof(PetscVTUReal) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0);
291: #else
292:             PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset));
293:             boffset += gpiece[r].ncells * 3 * sizeof(PetscVTUReal) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0);
294: #endif
295:             i += fbs;
296:           } else {
297:             for (j = 0; j < fbs; j++) {
298:               const char *compName = NULL;
299:               char        finalname[256];
300:               if (fv) PetscCall(PetscFVGetComponentName(fv, j, &compName));
301:               if (compName) {
302:                 PetscCall(PetscSNPrintf(finalname, 255, "%s%s.%s", vecname, fieldname, compName));
303:               } else if (fbs > 1) {
304:                 PetscCall(PetscSNPrintf(finalname, 255, "%s%s.%" PetscInt_FMT, vecname, fieldname, j));
305:               } else {
306:                 PetscCall(PetscSNPrintf(finalname, 255, "%s%s", vecname, fieldname));
307:               }
308: #if defined(PETSC_USE_COMPLEX)
309:               PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s.Re\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset));
310:               boffset += gpiece[r].ncells * sizeof(PetscVTUReal) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0);
311:               PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s.Im\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset));
312:               boffset += gpiece[r].ncells * sizeof(PetscVTUReal) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0);
313: #else
314:               PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset));
315:               boffset += gpiece[r].ncells * sizeof(PetscVTUReal) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0);
316: #endif
317:               i++;
318:             }
319:           }
320:         }
321:         //PetscCheck(i == bs,comm,PETSC_ERR_PLIB,"Total number of field components %" PetscInt_FMT " != block size %" PetscInt_FMT,i,bs);
322:       }
323:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "      </CellData>\n"));

325:       /*
326:        * Point Data headers
327:        */
328:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "      <PointData>\n"));
329:       for (link = vtk->link; link; link = link->next) {
330:         Vec          X = (Vec)link->vec;
331:         DM           dmX;
332:         PetscInt     bs      = 1, nfields, field;
333:         const char  *vecname = "";
334:         PetscSection section;
335:         if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
336:         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. */
337:           PetscCall(PetscObjectGetName((PetscObject)X, &vecname));
338:         }
339:         PetscCall(VecGetDM(X, &dmX));
340:         if (!dmX) dmX = dm;
341:         PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)&section));
342:         if (!section) PetscCall(DMGetLocalSection(dmX, &section));
343:         if (vEnd > vStart) PetscCall(PetscSectionGetDof(section, vStart, &bs));
344:         PetscCall(PetscSectionGetNumFields(section, &nfields));
345:         field = 0;
346:         if (link->field >= 0) {
347:           field   = link->field;
348:           nfields = field + 1;
349:         }
350:         for (i = 0; field < (nfields ? nfields : 1); field++) {
351:           PetscInt    fbs, j;
352:           const char *fieldname = NULL;
353:           char        buf[256];
354:           if (nfields) { /* We have user-defined fields/components */
355:             PetscCall(PetscSectionGetFieldDof(section, vStart, field, &fbs));
356:             PetscCall(PetscSectionGetFieldName(section, field, &fieldname));
357:           } else fbs = bs; /* Say we have one field with 'bs' components */
358:           if (nfields && !fieldname) {
359:             PetscCall(PetscSNPrintf(buf, sizeof(buf), "PointField%" PetscInt_FMT, field));
360:             fieldname = buf;
361:           }
362:           if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) {
363:             PetscCheck(fbs <= 3, PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_SIZ, "Point vector fields can have at most 3 components, %" PetscInt_FMT " given", fbs);
364: #if defined(PETSC_USE_COMPLEX)
365:             PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s%s.Re\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset));
366:             boffset += gpiece[r].nvertices * 3 * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0);
367:             PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s%s.Im\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset));
368:             boffset += gpiece[r].nvertices * 3 * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0);
369: #else
370:             PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset));
371:             boffset += gpiece[r].nvertices * 3 * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0);
372: #endif
373:           } else {
374:             for (j = 0; j < fbs; j++) {
375:               const char *compName = NULL;
376:               char        finalname[256];
377:               PetscCall(PetscSectionGetComponentName(section, field, j, &compName));
378:               PetscCall(PetscSNPrintf(finalname, 255, "%s%s.%s", vecname, fieldname, compName));
379: #if defined(PETSC_USE_COMPLEX)
380:               PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s.Re\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset));
381:               boffset += gpiece[r].nvertices * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0);
382:               PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s.Im\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset));
383:               boffset += gpiece[r].nvertices * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0);
384: #else
385:               PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset));
386:               boffset += gpiece[r].nvertices * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0);
387: #endif
388:             }
389:           }
390:         }
391:       }
392:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "      </PointData>\n"));
393:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "    </Piece>\n"));
394:     }
395:   }

397:   PetscCall(PetscFPrintf(comm, fp, "  </UnstructuredGrid>\n"));
398:   PetscCall(PetscFPrintf(comm, fp, "  <AppendedData encoding=\"raw\">\n"));
399:   PetscCall(PetscFPrintf(comm, fp, "_"));

401:   if (rank == 0) {
402:     PetscInt maxsize = 0;
403:     for (r = 0; r < size; r++) {
404:       maxsize = PetscMax(maxsize, (PetscInt)(gpiece[r].nvertices * 3 * sizeof(PetscVTUReal)));
405:       maxsize = PetscMax(maxsize, (PetscInt)(gpiece[r].ncells * 3 * sizeof(PetscVTUReal)));
406:       maxsize = PetscMax(maxsize, (PetscInt)(gpiece[r].nconn * sizeof(PetscVTKInt)));
407:     }
408:     PetscCall(PetscMalloc(maxsize, &buffer));
409:   }
410:   for (r = 0; r < size; r++) {
411:     if (r == rank) {
412:       PetscInt nsend;
413:       { /* Position */
414:         const PetscScalar *x, *cx = NULL;
415:         PetscVTUReal      *y = NULL;
416:         Vec                coords, cellCoords;
417:         PetscBool          copy;

419:         PetscCall(DMGetCoordinatesLocal(dm, &coords));
420:         PetscCall(VecGetArrayRead(coords, &x));
421:         PetscCall(DMGetCellCoordinatesLocal(dm, &cellCoords));
422:         if (cellCoords) PetscCall(VecGetArrayRead(cellCoords, &cx));
423: #if defined(PETSC_USE_COMPLEX)
424:         copy = PETSC_TRUE;
425: #else
426:         copy = (PetscBool)(dimEmbed != 3 || localized || (sizeof(PetscReal) != sizeof(PetscVTUReal)));
427: #endif
428:         if (copy) {
429:           PetscCall(PetscMalloc1(piece.nvertices * 3, &y));
430:           if (localized) {
431:             PetscInt cnt;
432:             for (c = cStart, cnt = 0; c < cEnd; c++) {
433:               PetscInt off, dof;

435:               PetscCall(PetscSectionGetDof(cellCoordSection, c, &dof));
436:               if (!dof) {
437:                 PetscInt *closure = NULL;
438:                 PetscInt  closureSize;

440:                 PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
441:                 for (v = 0; v < closureSize * 2; v += 2) {
442:                   if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
443:                     PetscCall(PetscSectionGetOffset(coordSection, closure[v], &off));
444:                     if (dimEmbed != 3) {
445:                       y[cnt * 3 + 0] = (PetscVTUReal)PetscRealPart(x[off + 0]);
446:                       y[cnt * 3 + 1] = (PetscVTUReal)((dimEmbed > 1) ? PetscRealPart(x[off + 1]) : 0.0);
447:                       y[cnt * 3 + 2] = (PetscVTUReal)0.0;
448:                     } else {
449:                       y[cnt * 3 + 0] = (PetscVTUReal)PetscRealPart(x[off + 0]);
450:                       y[cnt * 3 + 1] = (PetscVTUReal)PetscRealPart(x[off + 1]);
451:                       y[cnt * 3 + 2] = (PetscVTUReal)PetscRealPart(x[off + 2]);
452:                     }
453:                     cnt++;
454:                   }
455:                 }
456:                 PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
457:               } else {
458:                 PetscCall(PetscSectionGetOffset(cellCoordSection, c, &off));
459:                 if (dimEmbed != 3) {
460:                   for (i = 0; i < dof / dimEmbed; i++) {
461:                     y[cnt * 3 + 0] = (PetscVTUReal)PetscRealPart(cx[off + i * dimEmbed + 0]);
462:                     y[cnt * 3 + 1] = (PetscVTUReal)((dimEmbed > 1) ? PetscRealPart(cx[off + i * dimEmbed + 1]) : 0.0);
463:                     y[cnt * 3 + 2] = (PetscVTUReal)0.0;
464:                     cnt++;
465:                   }
466:                 } else {
467:                   for (i = 0; i < dof; i++) y[cnt * 3 + i] = (PetscVTUReal)PetscRealPart(cx[off + i]);
468:                   cnt += dof / dimEmbed;
469:                 }
470:               }
471:             }
472:             PetscCheck(cnt == piece.nvertices, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count does not match");
473:           } else {
474:             for (i = 0; i < piece.nvertices; i++) {
475:               y[i * 3 + 0] = (PetscVTUReal)PetscRealPart(x[i * dimEmbed + 0]);
476:               y[i * 3 + 1] = (PetscVTUReal)((dimEmbed > 1) ? PetscRealPart(x[i * dimEmbed + 1]) : 0.);
477:               y[i * 3 + 2] = (PetscVTUReal)((dimEmbed > 2) ? PetscRealPart(x[i * dimEmbed + 2]) : 0.);
478:             }
479:           }
480:         }
481:         nsend = piece.nvertices * 3;
482:         PetscCall(TransferWrite(comm, viewer, fp, r, 0, copy ? (const void *)y : (const void *)x, buffer, nsend, MPIU_VTUREAL, tag));
483:         PetscCall(PetscFree(y));
484:         PetscCall(VecRestoreArrayRead(coords, &x));
485:         if (cellCoords) PetscCall(VecRestoreArrayRead(cellCoords, &cx));
486:       }
487:       { /* Connectivity, offsets, types */
488:         PetscVTKInt  *connectivity = NULL, *offsets = NULL;
489:         PetscVTKType *types = NULL;
490:         PetscCall(DMPlexGetVTKConnectivity(dm, localized, &piece, &connectivity, &offsets, &types));
491:         PetscCall(TransferWrite(comm, viewer, fp, r, 0, connectivity, buffer, piece.nconn, MPI_INT, tag));
492:         PetscCall(TransferWrite(comm, viewer, fp, r, 0, offsets, buffer, piece.ncells, MPI_INT, tag));
493:         PetscCall(TransferWrite(comm, viewer, fp, r, 0, types, buffer, piece.ncells, MPI_CHAR, tag));
494:         PetscCall(PetscFree3(connectivity, offsets, types));
495:       }
496:       { /* Owners (cell data) */
497:         PetscVTKInt *owners;
498:         PetscCall(PetscMalloc1(piece.ncells, &owners));
499:         for (i = 0; i < piece.ncells; i++) owners[i] = rank;
500:         PetscCall(TransferWrite(comm, viewer, fp, r, 0, owners, buffer, piece.ncells, MPI_INT, tag));
501:         PetscCall(PetscFree(owners));
502:       }
503:       /* Cell data */
504:       for (link = vtk->link; link; link = link->next) {
505:         Vec                X = (Vec)link->vec;
506:         DM                 dmX;
507:         const PetscScalar *x;
508:         PetscVTUReal      *y;
509:         PetscInt           bs      = 1, nfields, field;
510:         PetscSection       section = NULL;

512:         if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
513:         PetscCall(VecGetDM(X, &dmX));
514:         if (!dmX) dmX = dm;
515:         PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)&section));
516:         if (!section) PetscCall(DMGetLocalSection(dmX, &section));
517:         if (cEnd > cStart) PetscCall(PetscSectionGetDof(section, cStart, &bs));
518:         PetscCall(PetscSectionGetNumFields(section, &nfields));
519:         field = 0;
520:         if (link->field >= 0) {
521:           field   = link->field;
522:           nfields = field + 1;
523:         }
524:         PetscCall(VecGetArrayRead(X, &x));
525:         PetscCall(PetscMalloc1(piece.ncells * 3, &y));
526:         for (i = 0; field < (nfields ? nfields : 1); field++) {
527:           PetscInt     fbs, j;
528:           PetscFV      fv = NULL;
529:           PetscObject  f;
530:           PetscClassId fClass;
531:           PetscBool    vector;
532:           if (nfields && cEnd > cStart) { /* We have user-defined fields/components */
533:             PetscCall(PetscSectionGetFieldDof(section, cStart, field, &fbs));
534:           } else fbs = bs; /* Say we have one field with 'bs' components */
535:           PetscCall(DMGetField(dmX, field, NULL, &f));
536:           PetscCall(PetscObjectGetClassId(f, &fClass));
537:           if (fClass == PETSCFV_CLASSID) fv = (PetscFV)f;
538:           vector = PETSC_FALSE;
539:           if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) {
540:             vector = PETSC_TRUE;
541:             for (j = 0; j < fbs; j++) {
542:               const char *compName = NULL;
543:               if (fv) {
544:                 PetscCall(PetscFVGetComponentName(fv, j, &compName));
545:                 if (compName) break;
546:               }
547:             }
548:             if (j < fbs) vector = PETSC_FALSE;
549:           }
550:           if (vector) {
551:             PetscInt cnt, l;
552:             for (l = 0; l < loops_per_scalar; l++) {
553:               for (c = cStart, cnt = 0; c < cEnd; c++) {
554:                 const PetscScalar *xpoint;
555:                 PetscInt           off, j;

557:                 if (hasLabel) { /* Ignore some cells */
558:                   PetscInt value;
559:                   PetscCall(DMGetLabelValue(dmX, "vtk", c, &value));
560:                   if (value != 1) continue;
561:                 }
562:                 if (nfields) {
563:                   PetscCall(PetscSectionGetFieldOffset(section, c, field, &off));
564:                 } else {
565:                   PetscCall(PetscSectionGetOffset(section, c, &off));
566:                 }
567:                 xpoint = &x[off];
568:                 for (j = 0; j < fbs; j++) y[cnt++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[j]) : PetscRealPart(xpoint[j]));
569:                 for (; j < 3; j++) y[cnt++] = 0.;
570:               }
571:               PetscCheck(cnt == piece.ncells * 3, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count does not match");
572:               PetscCall(TransferWrite(comm, viewer, fp, r, 0, y, buffer, piece.ncells * 3, MPIU_VTUREAL, tag));
573:             }
574:           } else {
575:             for (i = 0; i < fbs; i++) {
576:               PetscInt cnt, l;
577:               for (l = 0; l < loops_per_scalar; l++) {
578:                 for (c = cStart, cnt = 0; c < cEnd; c++) {
579:                   const PetscScalar *xpoint;
580:                   PetscInt           off;

582:                   if (hasLabel) { /* Ignore some cells */
583:                     PetscInt value;
584:                     PetscCall(DMGetLabelValue(dmX, "vtk", c, &value));
585:                     if (value != 1) continue;
586:                   }
587:                   if (nfields) {
588:                     PetscCall(PetscSectionGetFieldOffset(section, c, field, &off));
589:                   } else {
590:                     PetscCall(PetscSectionGetOffset(section, c, &off));
591:                   }
592:                   xpoint   = &x[off];
593:                   y[cnt++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i]));
594:                 }
595:                 PetscCheck(cnt == piece.ncells, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count does not match");
596:                 PetscCall(TransferWrite(comm, viewer, fp, r, 0, y, buffer, piece.ncells, MPIU_VTUREAL, tag));
597:               }
598:             }
599:           }
600:         }
601:         PetscCall(PetscFree(y));
602:         PetscCall(VecRestoreArrayRead(X, &x));
603:       }
604:       /* point data */
605:       for (link = vtk->link; link; link = link->next) {
606:         Vec                X = (Vec)link->vec;
607:         DM                 dmX;
608:         const PetscScalar *x;
609:         PetscVTUReal      *y;
610:         PetscInt           bs      = 1, nfields, field;
611:         PetscSection       section = NULL;

613:         if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
614:         PetscCall(VecGetDM(X, &dmX));
615:         if (!dmX) dmX = dm;
616:         PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)&section));
617:         if (!section) PetscCall(DMGetLocalSection(dmX, &section));
618:         if (vEnd > vStart) PetscCall(PetscSectionGetDof(section, vStart, &bs));
619:         PetscCall(PetscSectionGetNumFields(section, &nfields));
620:         field = 0;
621:         if (link->field >= 0) {
622:           field   = link->field;
623:           nfields = field + 1;
624:         }
625:         PetscCall(VecGetArrayRead(X, &x));
626:         PetscCall(PetscMalloc1(piece.nvertices * 3, &y));
627:         for (i = 0; field < (nfields ? nfields : 1); field++) {
628:           PetscInt fbs, j;
629:           if (nfields && vEnd > vStart) { /* We have user-defined fields/components */
630:             PetscCall(PetscSectionGetFieldDof(section, vStart, field, &fbs));
631:           } else fbs = bs; /* Say we have one field with 'bs' components */
632:           if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) {
633:             PetscInt cnt, l;
634:             for (l = 0; l < loops_per_scalar; l++) {
635:               if (!localized) {
636:                 for (v = vStart, cnt = 0; v < vEnd; v++) {
637:                   PetscInt           off;
638:                   const PetscScalar *xpoint;

640:                   if (nfields) {
641:                     PetscCall(PetscSectionGetFieldOffset(section, v, field, &off));
642:                   } else {
643:                     PetscCall(PetscSectionGetOffset(section, v, &off));
644:                   }
645:                   xpoint = &x[off];
646:                   for (j = 0; j < fbs; j++) y[cnt++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[j]) : PetscRealPart(xpoint[j]));
647:                   for (; j < 3; j++) y[cnt++] = 0.;
648:                 }
649:               } else {
650:                 for (c = cStart, cnt = 0; c < cEnd; c++) {
651:                   PetscInt *closure = NULL;
652:                   PetscInt  closureSize, off;

654:                   PetscCall(DMPlexGetTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure));
655:                   for (v = 0, off = 0; v < closureSize * 2; v += 2) {
656:                     if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
657:                       PetscInt           voff;
658:                       const PetscScalar *xpoint;

660:                       if (nfields) {
661:                         PetscCall(PetscSectionGetFieldOffset(section, closure[v], field, &voff));
662:                       } else {
663:                         PetscCall(PetscSectionGetOffset(section, closure[v], &voff));
664:                       }
665:                       xpoint = &x[voff];
666:                       for (j = 0; j < fbs; j++) y[cnt + off++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[j]) : PetscRealPart(xpoint[j]));
667:                       for (; j < 3; j++) y[cnt + off++] = 0.;
668:                     }
669:                   }
670:                   cnt += off;
671:                   PetscCall(DMPlexRestoreTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure));
672:                 }
673:               }
674:               PetscCheck(cnt == piece.nvertices * 3, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count does not match");
675:               PetscCall(TransferWrite(comm, viewer, fp, r, 0, y, buffer, piece.nvertices * 3, MPIU_VTUREAL, tag));
676:             }
677:           } else {
678:             for (i = 0; i < fbs; i++) {
679:               PetscInt cnt, l;
680:               for (l = 0; l < loops_per_scalar; l++) {
681:                 if (!localized) {
682:                   for (v = vStart, cnt = 0; v < vEnd; v++) {
683:                     PetscInt           off;
684:                     const PetscScalar *xpoint;

686:                     if (nfields) {
687:                       PetscCall(PetscSectionGetFieldOffset(section, v, field, &off));
688:                     } else {
689:                       PetscCall(PetscSectionGetOffset(section, v, &off));
690:                     }
691:                     xpoint   = &x[off];
692:                     y[cnt++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i]));
693:                   }
694:                 } else {
695:                   for (c = cStart, cnt = 0; c < cEnd; c++) {
696:                     PetscInt *closure = NULL;
697:                     PetscInt  closureSize, off;

699:                     PetscCall(DMPlexGetTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure));
700:                     for (v = 0, off = 0; v < closureSize * 2; v += 2) {
701:                       if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
702:                         PetscInt           voff;
703:                         const PetscScalar *xpoint;

705:                         if (nfields) {
706:                           PetscCall(PetscSectionGetFieldOffset(section, closure[v], field, &voff));
707:                         } else {
708:                           PetscCall(PetscSectionGetOffset(section, closure[v], &voff));
709:                         }
710:                         xpoint         = &x[voff];
711:                         y[cnt + off++] = (l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i]));
712:                       }
713:                     }
714:                     cnt += off;
715:                     PetscCall(DMPlexRestoreTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure));
716:                   }
717:                 }
718:                 PetscCheck(cnt == piece.nvertices, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count does not match");
719:                 PetscCall(TransferWrite(comm, viewer, fp, r, 0, y, buffer, piece.nvertices, MPIU_VTUREAL, tag));
720:               }
721:             }
722:           }
723:         }
724:         PetscCall(PetscFree(y));
725:         PetscCall(VecRestoreArrayRead(X, &x));
726:       }
727:     } else if (rank == 0) {
728:       PetscInt l;

730:       PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].nvertices * 3, MPIU_VTUREAL, tag)); /* positions */
731:       PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].nconn, MPI_INT, tag));              /* connectivity */
732:       PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells, MPI_INT, tag));             /* offsets */
733:       PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells, MPI_CHAR, tag));            /* types */
734:       PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells, MPI_INT, tag));             /* owner rank (cells) */
735:       /* all cell data */
736:       for (link = vtk->link; link; link = link->next) {
737:         Vec          X  = (Vec)link->vec;
738:         PetscInt     bs = 1, nfields, field;
739:         DM           dmX;
740:         PetscSection section = NULL;

742:         if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
743:         PetscCall(VecGetDM(X, &dmX));
744:         if (!dmX) dmX = dm;
745:         PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)&section));
746:         if (!section) PetscCall(DMGetLocalSection(dmX, &section));
747:         if (cEnd > cStart) PetscCall(PetscSectionGetDof(section, cStart, &bs));
748:         PetscCall(PetscSectionGetNumFields(section, &nfields));
749:         field = 0;
750:         if (link->field >= 0) {
751:           field   = link->field;
752:           nfields = field + 1;
753:         }
754:         for (i = 0; field < (nfields ? nfields : 1); field++) {
755:           PetscInt     fbs, j;
756:           PetscFV      fv = NULL;
757:           PetscObject  f;
758:           PetscClassId fClass;
759:           PetscBool    vector;
760:           if (nfields && cEnd > cStart) { /* We have user-defined fields/components */
761:             PetscCall(PetscSectionGetFieldDof(section, cStart, field, &fbs));
762:           } else fbs = bs; /* Say we have one field with 'bs' components */
763:           PetscCall(DMGetField(dmX, field, NULL, &f));
764:           PetscCall(PetscObjectGetClassId(f, &fClass));
765:           if (fClass == PETSCFV_CLASSID) fv = (PetscFV)f;
766:           vector = PETSC_FALSE;
767:           if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) {
768:             vector = PETSC_TRUE;
769:             for (j = 0; j < fbs; j++) {
770:               const char *compName = NULL;
771:               if (fv) {
772:                 PetscCall(PetscFVGetComponentName(fv, j, &compName));
773:                 if (compName) break;
774:               }
775:             }
776:             if (j < fbs) vector = PETSC_FALSE;
777:           }
778:           if (vector) {
779:             for (l = 0; l < loops_per_scalar; l++) PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells * 3, MPIU_VTUREAL, tag));
780:           } else {
781:             for (i = 0; i < fbs; i++) {
782:               for (l = 0; l < loops_per_scalar; l++) PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells, MPIU_VTUREAL, tag));
783:             }
784:           }
785:         }
786:       }
787:       /* all point data */
788:       for (link = vtk->link; link; link = link->next) {
789:         Vec          X = (Vec)link->vec;
790:         DM           dmX;
791:         PetscInt     bs      = 1, nfields, field;
792:         PetscSection section = NULL;

794:         if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
795:         PetscCall(VecGetDM(X, &dmX));
796:         if (!dmX) dmX = dm;
797:         PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)&section));
798:         if (!section) PetscCall(DMGetLocalSection(dmX, &section));
799:         if (vEnd > vStart) PetscCall(PetscSectionGetDof(section, vStart, &bs));
800:         PetscCall(PetscSectionGetNumFields(section, &nfields));
801:         field = 0;
802:         if (link->field >= 0) {
803:           field   = link->field;
804:           nfields = field + 1;
805:         }
806:         for (i = 0; field < (nfields ? nfields : 1); field++) {
807:           PetscInt fbs;
808:           if (nfields && vEnd > vStart) { /* We have user-defined fields/components */
809:             PetscCall(PetscSectionGetFieldDof(section, vStart, field, &fbs));
810:           } else fbs = bs; /* Say we have one field with 'bs' components */
811:           if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) {
812:             for (l = 0; l < loops_per_scalar; l++) PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].nvertices * 3, MPIU_VTUREAL, tag));
813:           } else {
814:             for (i = 0; i < fbs; i++) {
815:               for (l = 0; l < loops_per_scalar; l++) PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].nvertices, MPIU_VTUREAL, tag));
816:             }
817:           }
818:         }
819:       }
820:     }
821:   }
822:   PetscCall(PetscFree(gpiece));
823:   PetscCall(PetscFree(buffer));
824:   PetscCall(PetscFPrintf(comm, fp, "\n  </AppendedData>\n"));
825:   PetscCall(PetscFPrintf(comm, fp, "</VTKFile>\n"));
826:   PetscCall(PetscFClose(comm, fp));
827:   PetscFunctionReturn(PETSC_SUCCESS);
828: }