1 #include <petsc/private/dmpleximpl.h>
2 #include <../src/sys/classes/viewer/impls/vtk/vtkvimpl.h>
3
4 typedef struct {
5 PetscInt nvertices;
6 PetscInt ncells;
7 PetscInt nconn; /* number of entries in cell->vertex connectivity array */
8 } PieceInfo;
9
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
25
TransferWrite(MPI_Comm comm,PetscViewer viewer,FILE * fp,PetscMPIInt srank,PetscMPIInt root,const void * send,void * recv,PetscCount count,MPI_Datatype mpidatatype,PetscMPIInt tag)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;
29
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 }
50
DMPlexGetVTKConnectivity(DM dm,PetscBool localized,PieceInfo * piece,PetscVTKInt ** oconn,PetscVTKInt ** ooffsets,PetscVTKType ** otypes)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;
57
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;
69
70 countcell = 0;
71 countconn = 0;
72 for (c = cStart; c < cEnd; ++c) {
73 PetscInt nverts, dof = 0, celltype, startoffset, nC = 0;
74
75 if (hasLabel) {
76 PetscInt value;
77
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;
86
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 }
99
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]));
107
108 nverts = countconn - startoffset;
109 PetscCall(DMPlexVTKGetCellType_Internal(dm, dim, nverts, &celltype));
110
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 }
121
DMPlexGetNonEmptyComm_Private(DM dm,MPI_Comm * comm)122 PETSC_INTERN PetscErrorCode DMPlexGetNonEmptyComm_Private(DM dm, MPI_Comm *comm)
123 {
124 DM_Plex *mesh = (DM_Plex *)dm->data;
125
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;
131
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 }
145
DMGetFieldIfFV_Private(DM dm,PetscInt field,PetscFV * fv)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;
151
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 }
162
163 /*
164 Write all fields that have been provided to the viewer
165 Multi-block XML format with binary appended data.
166 */
DMPlexVTKWriteAll_VTU(DM dm,PetscViewer viewer)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;
181
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));
192
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));
197
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"));
202
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;
212
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;
220
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));
237
238 /*
239 * Write file header
240 */
241 if (rank == 0) {
242 PetscInt64 boffset = 0;
243
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"));
260
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;
295
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"));
357
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 }
429
430 PetscCall(PetscFPrintf(comm, fp, " </UnstructuredGrid>\n"));
431 PetscCall(PetscFPrintf(comm, fp, " <AppendedData encoding=\"raw\">\n"));
432 PetscCall(PetscFPrintf(comm, fp, "_"));
433
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)));
451
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;
462
463 PetscCall(PetscSectionGetDof(cellCoordSection, c, &dof));
464 if (!dof) {
465 PetscInt *closure = NULL;
466 PetscInt closureSize;
467
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;
527
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;
542
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;
561
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;
584
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;
609
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;
640
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;
667
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;
681
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;
687
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;
713
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;
726
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;
732
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;
757
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;
769
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;
786
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;
818
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 }
861