Actual source code: plexvtk.c

  1: #define PETSCDM_DLL
  2: #include <petsc/private/dmpleximpl.h>
  3: #include <../src/sys/classes/viewer/impls/vtk/vtkvimpl.h>

  5: PetscErrorCode DMPlexVTKGetCellType_Internal(DM dm, PetscInt dim, PetscInt corners, PetscInt *cellType)
  6: {
  7:   *cellType = -1;
  8:   switch (dim) {
  9:   case 0:
 10:     switch (corners) {
 11:     case 1:
 12:       *cellType = 1; /* VTK_VERTEX */
 13:       break;
 14:     default:
 15:       break;
 16:     }
 17:     break;
 18:   case 1:
 19:     switch (corners) {
 20:     case 2:
 21:       *cellType = 3; /* VTK_LINE */
 22:       break;
 23:     case 3:
 24:       *cellType = 21; /* VTK_QUADRATIC_EDGE */
 25:       break;
 26:     default:
 27:       break;
 28:     }
 29:     break;
 30:   case 2:
 31:     switch (corners) {
 32:     case 3:
 33:       *cellType = 5; /* VTK_TRIANGLE */
 34:       break;
 35:     case 4:
 36:       *cellType = 9; /* VTK_QUAD */
 37:       break;
 38:     case 6:
 39:       *cellType = 22; /* VTK_QUADRATIC_TRIANGLE */
 40:       break;
 41:     case 9:
 42:       *cellType = 23; /* VTK_QUADRATIC_QUAD */
 43:       break;
 44:     default:
 45:       break;
 46:     }
 47:     break;
 48:   case 3:
 49:     switch (corners) {
 50:     case 4:
 51:       *cellType = 10; /* VTK_TETRA */
 52:       break;
 53:     case 5:
 54:       *cellType = 14; /* VTK_PYRAMID */
 55:       break;
 56:     case 6:
 57:       *cellType = 13; /* VTK_WEDGE */
 58:       break;
 59:     case 8:
 60:       *cellType = 12; /* VTK_HEXAHEDRON */
 61:       break;
 62:     case 10:
 63:       *cellType = 24; /* VTK_QUADRATIC_TETRA */
 64:       break;
 65:     case 27:
 66:       *cellType = 29; /* VTK_QUADRATIC_HEXAHEDRON */
 67:       break;
 68:     default:
 69:       break;
 70:     }
 71:   }
 72:   return 0;
 73: }

 75: static PetscErrorCode DMPlexVTKWriteCells_ASCII(DM dm, FILE *fp, PetscInt *totalCells)
 76: {
 77:   MPI_Comm        comm;
 78:   DMLabel         label;
 79:   IS              globalVertexNumbers = NULL;
 80:   const PetscInt *gvertex;
 81:   PetscInt        dim;
 82:   PetscInt        numCorners = 0, totCorners = 0, maxCorners, *corners;
 83:   PetscInt        numCells = 0, totCells = 0, maxCells, cellHeight;
 84:   PetscInt        numLabelCells, maxLabelCells, cStart, cEnd, c, vStart, vEnd, v;
 85:   PetscMPIInt     size, rank, proc, tag;

 87:   PetscObjectGetComm((PetscObject)dm, &comm);
 88:   PetscCommGetNewTag(comm, &tag);
 89:   MPI_Comm_size(comm, &size);
 90:   MPI_Comm_rank(comm, &rank);
 91:   DMGetDimension(dm, &dim);
 92:   DMPlexGetVTKCellHeight(dm, &cellHeight);
 93:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
 94:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
 95:   DMGetLabel(dm, "vtk", &label);
 96:   DMGetStratumSize(dm, "vtk", 1, &numLabelCells);
 97:   MPIU_Allreduce(&numLabelCells, &maxLabelCells, 1, MPIU_INT, MPI_MAX, comm);
 98:   if (!maxLabelCells) label = NULL;
 99:   for (c = cStart; c < cEnd; ++c) {
100:     PetscInt *closure = NULL;
101:     PetscInt  closureSize, value;

103:     if (label) {
104:       DMLabelGetValue(label, c, &value);
105:       if (value != 1) continue;
106:     }
107:     DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
108:     for (v = 0; v < closureSize * 2; v += 2) {
109:       if ((closure[v] >= vStart) && (closure[v] < vEnd)) ++numCorners;
110:     }
111:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
112:     ++numCells;
113:   }
114:   maxCells = numCells;
115:   MPI_Reduce(&numCells, &totCells, 1, MPIU_INT, MPI_SUM, 0, comm);
116:   MPI_Reduce(&numCells, &maxCells, 1, MPIU_INT, MPI_MAX, 0, comm);
117:   MPI_Reduce(&numCorners, &totCorners, 1, MPIU_INT, MPI_SUM, 0, comm);
118:   MPI_Reduce(&numCorners, &maxCorners, 1, MPIU_INT, MPI_MAX, 0, comm);
119:   DMPlexGetVertexNumbering(dm, &globalVertexNumbers);
120:   ISGetIndices(globalVertexNumbers, &gvertex);
121:   PetscMalloc1(maxCells, &corners);
122:   PetscFPrintf(comm, fp, "CELLS %" PetscInt_FMT " %" PetscInt_FMT "\n", totCells, totCorners + totCells);
123:   if (rank == 0) {
124:     PetscInt *remoteVertices, *vertices;

126:     PetscMalloc1(maxCorners, &vertices);
127:     for (c = cStart, numCells = 0; c < cEnd; ++c) {
128:       PetscInt *closure = NULL;
129:       PetscInt  closureSize, value, nC = 0;

131:       if (label) {
132:         DMLabelGetValue(label, c, &value);
133:         if (value != 1) continue;
134:       }
135:       DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
136:       for (v = 0; v < closureSize * 2; v += 2) {
137:         if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
138:           const PetscInt gv = gvertex[closure[v] - vStart];
139:           vertices[nC++]    = gv < 0 ? -(gv + 1) : gv;
140:         }
141:       }
142:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
143:       DMPlexReorderCell(dm, c, vertices);
144:       corners[numCells++] = nC;
145:       PetscFPrintf(comm, fp, "%" PetscInt_FMT " ", nC);
146:       for (v = 0; v < nC; ++v) PetscFPrintf(comm, fp, " %" PetscInt_FMT, vertices[v]);
147:       PetscFPrintf(comm, fp, "\n");
148:     }
149:     if (size > 1) PetscMalloc1(maxCorners + maxCells, &remoteVertices);
150:     for (proc = 1; proc < size; ++proc) {
151:       MPI_Status status;

153:       MPI_Recv(&numCorners, 1, MPIU_INT, proc, tag, comm, &status);
154:       MPI_Recv(remoteVertices, numCorners, MPIU_INT, proc, tag, comm, &status);
155:       for (c = 0; c < numCorners;) {
156:         PetscInt nC = remoteVertices[c++];

158:         for (v = 0; v < nC; ++v, ++c) vertices[v] = remoteVertices[c];
159:         PetscFPrintf(comm, fp, "%" PetscInt_FMT " ", nC);
160:         for (v = 0; v < nC; ++v) PetscFPrintf(comm, fp, " %" PetscInt_FMT, vertices[v]);
161:         PetscFPrintf(comm, fp, "\n");
162:       }
163:     }
164:     if (size > 1) PetscFree(remoteVertices);
165:     PetscFree(vertices);
166:   } else {
167:     PetscInt *localVertices, numSend = numCells + numCorners, k = 0;

169:     PetscMalloc1(numSend, &localVertices);
170:     for (c = cStart, numCells = 0; c < cEnd; ++c) {
171:       PetscInt *closure = NULL;
172:       PetscInt  closureSize, value, nC = 0;

174:       if (label) {
175:         DMLabelGetValue(label, c, &value);
176:         if (value != 1) continue;
177:       }
178:       DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
179:       for (v = 0; v < closureSize * 2; v += 2) {
180:         if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
181:           const PetscInt gv = gvertex[closure[v] - vStart];
182:           closure[nC++]     = gv < 0 ? -(gv + 1) : gv;
183:         }
184:       }
185:       corners[numCells++] = nC;
186:       localVertices[k++]  = nC;
187:       for (v = 0; v < nC; ++v, ++k) localVertices[k] = closure[v];
188:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
189:       DMPlexReorderCell(dm, c, localVertices + k - nC);
190:     }
192:     MPI_Send(&numSend, 1, MPIU_INT, 0, tag, comm);
193:     MPI_Send(localVertices, numSend, MPIU_INT, 0, tag, comm);
194:     PetscFree(localVertices);
195:   }
196:   ISRestoreIndices(globalVertexNumbers, &gvertex);
197:   PetscFPrintf(comm, fp, "CELL_TYPES %" PetscInt_FMT "\n", totCells);
198:   if (rank == 0) {
199:     PetscInt cellType;

201:     for (c = 0; c < numCells; ++c) {
202:       DMPlexVTKGetCellType_Internal(dm, dim, corners[c], &cellType);
203:       PetscFPrintf(comm, fp, "%" PetscInt_FMT "\n", cellType);
204:     }
205:     for (proc = 1; proc < size; ++proc) {
206:       MPI_Status status;

208:       MPI_Recv(&numCells, 1, MPIU_INT, proc, tag, comm, &status);
209:       MPI_Recv(corners, numCells, MPIU_INT, proc, tag, comm, &status);
210:       for (c = 0; c < numCells; ++c) {
211:         DMPlexVTKGetCellType_Internal(dm, dim, corners[c], &cellType);
212:         PetscFPrintf(comm, fp, "%" PetscInt_FMT "\n", cellType);
213:       }
214:     }
215:   } else {
216:     MPI_Send(&numCells, 1, MPIU_INT, 0, tag, comm);
217:     MPI_Send(corners, numCells, MPIU_INT, 0, tag, comm);
218:   }
219:   PetscFree(corners);
220:   *totalCells = totCells;
221:   return 0;
222: }

224: static PetscErrorCode DMPlexVTKWritePartition_ASCII(DM dm, FILE *fp)
225: {
226:   MPI_Comm    comm;
227:   PetscInt    numCells = 0, cellHeight;
228:   PetscInt    numLabelCells, cStart, cEnd, c;
229:   PetscMPIInt size, rank, proc, tag;
230:   PetscBool   hasLabel;

232:   PetscObjectGetComm((PetscObject)dm, &comm);
233:   PetscCommGetNewTag(comm, &tag);
234:   MPI_Comm_size(comm, &size);
235:   MPI_Comm_rank(comm, &rank);
236:   DMPlexGetVTKCellHeight(dm, &cellHeight);
237:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
238:   DMGetStratumSize(dm, "vtk", 1, &numLabelCells);
239:   hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;
240:   for (c = cStart; c < cEnd; ++c) {
241:     if (hasLabel) {
242:       PetscInt value;

244:       DMGetLabelValue(dm, "vtk", c, &value);
245:       if (value != 1) continue;
246:     }
247:     ++numCells;
248:   }
249:   if (rank == 0) {
250:     for (c = 0; c < numCells; ++c) PetscFPrintf(comm, fp, "%d\n", rank);
251:     for (proc = 1; proc < size; ++proc) {
252:       MPI_Status status;

254:       MPI_Recv(&numCells, 1, MPIU_INT, proc, tag, comm, &status);
255:       for (c = 0; c < numCells; ++c) PetscFPrintf(comm, fp, "%d\n", proc);
256:     }
257:   } else {
258:     MPI_Send(&numCells, 1, MPIU_INT, 0, tag, comm);
259:   }
260:   return 0;
261: }

263: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
264: typedef double PetscVTKReal;
265: #elif defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL___FP16)
266: typedef float PetscVTKReal;
267: #else
268: typedef PetscReal PetscVTKReal;
269: #endif

271: static PetscErrorCode DMPlexVTKWriteSection_ASCII(DM dm, PetscSection section, PetscSection globalSection, Vec v, FILE *fp, PetscInt enforceDof, PetscInt precision, PetscReal scale, PetscInt imag)
272: {
273:   MPI_Comm           comm;
274:   const MPI_Datatype mpiType = MPIU_SCALAR;
275:   PetscScalar       *array;
276:   PetscInt           numDof = 0, maxDof;
277:   PetscInt           numLabelCells, cellHeight, cStart, cEnd, numLabelVertices, vStart, vEnd, pStart, pEnd, p;
278:   PetscMPIInt        size, rank, proc, tag;
279:   PetscBool          hasLabel;

281:   PetscObjectGetComm((PetscObject)dm, &comm);
284:   if (precision < 0) precision = 6;
285:   PetscCommGetNewTag(comm, &tag);
286:   MPI_Comm_size(comm, &size);
287:   MPI_Comm_rank(comm, &rank);
288:   PetscSectionGetChart(section, &pStart, &pEnd);
289:   /* VTK only wants the values at cells or vertices */
290:   DMPlexGetVTKCellHeight(dm, &cellHeight);
291:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
292:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
293:   pStart = PetscMax(PetscMin(cStart, vStart), pStart);
294:   pEnd   = PetscMin(PetscMax(cEnd, vEnd), pEnd);
295:   DMGetStratumSize(dm, "vtk", 1, &numLabelCells);
296:   DMGetStratumSize(dm, "vtk", 2, &numLabelVertices);
297:   hasLabel = numLabelCells > 0 || numLabelVertices > 0 ? PETSC_TRUE : PETSC_FALSE;
298:   for (p = pStart; p < pEnd; ++p) {
299:     /* Reject points not either cells or vertices */
300:     if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
301:     if (hasLabel) {
302:       PetscInt value;

304:       if (((p >= cStart) && (p < cEnd) && numLabelCells) || ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
305:         DMGetLabelValue(dm, "vtk", p, &value);
306:         if (value != 1) continue;
307:       }
308:     }
309:     PetscSectionGetDof(section, p, &numDof);
310:     if (numDof) break;
311:   }
312:   MPIU_Allreduce(&numDof, &maxDof, 1, MPIU_INT, MPI_MAX, comm);
313:   enforceDof = PetscMax(enforceDof, maxDof);
314:   VecGetArray(v, &array);
315:   if (rank == 0) {
316:     PetscVTKReal dval;
317:     PetscScalar  val;
318:     char         formatString[8];

320:     PetscSNPrintf(formatString, 8, "%%.%" PetscInt_FMT "e", precision);
321:     for (p = pStart; p < pEnd; ++p) {
322:       /* Here we lose a way to filter points by keeping them out of the Numbering */
323:       PetscInt dof, off, goff, d;

325:       /* Reject points not either cells or vertices */
326:       if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
327:       if (hasLabel) {
328:         PetscInt value;

330:         if (((p >= cStart) && (p < cEnd) && numLabelCells) || ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
331:           DMGetLabelValue(dm, "vtk", p, &value);
332:           if (value != 1) continue;
333:         }
334:       }
335:       PetscSectionGetDof(section, p, &dof);
336:       PetscSectionGetOffset(section, p, &off);
337:       PetscSectionGetOffset(globalSection, p, &goff);
338:       if (dof && goff >= 0) {
339:         for (d = 0; d < dof; d++) {
340:           if (d > 0) PetscFPrintf(comm, fp, " ");
341:           val  = array[off + d];
342:           dval = (PetscVTKReal)((imag ? PetscImaginaryPart(val) : PetscRealPart(val)) * scale);
343:           PetscFPrintf(comm, fp, formatString, dval);
344:         }
345:         for (d = dof; d < enforceDof; d++) PetscFPrintf(comm, fp, " 0.0");
346:         PetscFPrintf(comm, fp, "\n");
347:       }
348:     }
349:     for (proc = 1; proc < size; ++proc) {
350:       PetscScalar *remoteValues;
351:       PetscInt     size = 0, d;
352:       MPI_Status   status;

354:       MPI_Recv(&size, 1, MPIU_INT, proc, tag, comm, &status);
355:       PetscMalloc1(size, &remoteValues);
356:       MPI_Recv(remoteValues, size, mpiType, proc, tag, comm, &status);
357:       for (p = 0; p < size / maxDof; ++p) {
358:         for (d = 0; d < maxDof; ++d) {
359:           if (d > 0) PetscFPrintf(comm, fp, " ");
360:           val  = remoteValues[p * maxDof + d];
361:           dval = (PetscVTKReal)((imag ? PetscImaginaryPart(val) : PetscRealPart(val)) * scale);
362:           PetscFPrintf(comm, fp, formatString, dval);
363:         }
364:         for (d = maxDof; d < enforceDof; ++d) PetscFPrintf(comm, fp, " 0.0");
365:         PetscFPrintf(comm, fp, "\n");
366:       }
367:       PetscFree(remoteValues);
368:     }
369:   } else {
370:     PetscScalar *localValues;
371:     PetscInt     size, k = 0;

373:     PetscSectionGetStorageSize(section, &size);
374:     PetscMalloc1(size, &localValues);
375:     for (p = pStart; p < pEnd; ++p) {
376:       PetscInt dof, off, goff, d;

378:       /* Reject points not either cells or vertices */
379:       if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
380:       if (hasLabel) {
381:         PetscInt value;

383:         if (((p >= cStart) && (p < cEnd) && numLabelCells) || ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
384:           DMGetLabelValue(dm, "vtk", p, &value);
385:           if (value != 1) continue;
386:         }
387:       }
388:       PetscSectionGetDof(section, p, &dof);
389:       PetscSectionGetOffset(section, p, &off);
390:       PetscSectionGetOffset(globalSection, p, &goff);
391:       if (goff >= 0) {
392:         for (d = 0; d < dof; ++d) localValues[k++] = array[off + d];
393:       }
394:     }
395:     MPI_Send(&k, 1, MPIU_INT, 0, tag, comm);
396:     MPI_Send(localValues, k, mpiType, 0, tag, comm);
397:     PetscFree(localValues);
398:   }
399:   VecRestoreArray(v, &array);
400:   return 0;
401: }

403: static PetscErrorCode DMPlexVTKWriteField_ASCII(DM dm, PetscSection section, PetscSection globalSection, Vec field, const char name[], FILE *fp, PetscInt enforceDof, PetscInt precision, PetscReal scale, PetscBool nameComplex, PetscInt imag)
404: {
405:   MPI_Comm comm;
406:   PetscInt numDof = 0, maxDof;
407:   PetscInt pStart, pEnd, p;

409:   PetscObjectGetComm((PetscObject)dm, &comm);
410:   PetscSectionGetChart(section, &pStart, &pEnd);
411:   for (p = pStart; p < pEnd; ++p) {
412:     PetscSectionGetDof(section, p, &numDof);
413:     if (numDof) break;
414:   }
415:   numDof = PetscMax(numDof, enforceDof);
416:   MPIU_Allreduce(&numDof, &maxDof, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
417:   if (!name) name = "Unknown";
418:   if (maxDof == 3) {
419:     if (nameComplex) {
420:       PetscFPrintf(comm, fp, "VECTORS %s.%s double\n", name, imag ? "Im" : "Re");
421:     } else {
422:       PetscFPrintf(comm, fp, "VECTORS %s double\n", name);
423:     }
424:   } else {
425:     if (nameComplex) {
426:       PetscFPrintf(comm, fp, "SCALARS %s.%s double %" PetscInt_FMT "\n", name, imag ? "Im" : "Re", maxDof);
427:     } else {
428:       PetscFPrintf(comm, fp, "SCALARS %s double %" PetscInt_FMT "\n", name, maxDof);
429:     }
430:     PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n");
431:   }
432:   DMPlexVTKWriteSection_ASCII(dm, section, globalSection, field, fp, enforceDof, precision, scale, imag);
433:   return 0;
434: }

436: static PetscErrorCode DMPlexVTKWriteAll_ASCII(DM dm, PetscViewer viewer)
437: {
438:   MPI_Comm                 comm;
439:   PetscViewer_VTK         *vtk = (PetscViewer_VTK *)viewer->data;
440:   FILE                    *fp;
441:   PetscViewerVTKObjectLink link;
442:   PetscInt                 totVertices, totCells = 0, loops_per_scalar, l;
443:   PetscBool                hasPoint = PETSC_FALSE, hasCell = PETSC_FALSE, writePartition = PETSC_FALSE, localized, writeComplex;
444:   const char              *dmname;

446: #if defined(PETSC_USE_COMPLEX)
447:   loops_per_scalar = 2;
448:   writeComplex     = PETSC_TRUE;
449: #else
450:   loops_per_scalar = 1;
451:   writeComplex     = PETSC_FALSE;
452: #endif
453:   DMGetCoordinatesLocalized(dm, &localized);
454:   PetscObjectGetComm((PetscObject)dm, &comm);
456:   PetscFOpen(comm, vtk->filename, "wb", &fp);
457:   PetscObjectGetName((PetscObject)dm, &dmname);
458:   PetscFPrintf(comm, fp, "# vtk DataFile Version 2.0\n");
459:   PetscFPrintf(comm, fp, "%s\n", dmname);
460:   PetscFPrintf(comm, fp, "ASCII\n");
461:   PetscFPrintf(comm, fp, "DATASET UNSTRUCTURED_GRID\n");
462:   /* Vertices */
463:   {
464:     PetscSection section, coordSection, globalCoordSection;
465:     Vec          coordinates;
466:     PetscReal    lengthScale;
467:     DMLabel      label;
468:     IS           vStratumIS;
469:     PetscLayout  vLayout;

471:     DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
472:     DMGetCoordinatesLocal(dm, &coordinates);
473:     DMPlexGetDepthLabel(dm, &label);
474:     DMLabelGetStratumIS(label, 0, &vStratumIS);
475:     DMGetCoordinateSection(dm, &section);                                   /* This section includes all points */
476:     PetscSectionCreateSubdomainSection(section, vStratumIS, &coordSection); /* This one includes just vertices */
477:     PetscSectionCreateGlobalSection(coordSection, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalCoordSection);
478:     PetscSectionGetPointLayout(comm, globalCoordSection, &vLayout);
479:     PetscLayoutGetSize(vLayout, &totVertices);
480:     PetscFPrintf(comm, fp, "POINTS %" PetscInt_FMT " double\n", totVertices);
481:     DMPlexVTKWriteSection_ASCII(dm, coordSection, globalCoordSection, coordinates, fp, 3, PETSC_DETERMINE, lengthScale, 0);
482:     ISDestroy(&vStratumIS);
483:     PetscLayoutDestroy(&vLayout);
484:     PetscSectionDestroy(&coordSection);
485:     PetscSectionDestroy(&globalCoordSection);
486:   }
487:   /* Cells */
488:   DMPlexVTKWriteCells_ASCII(dm, fp, &totCells);
489:   /* Vertex fields */
490:   for (link = vtk->link; link; link = link->next) {
491:     if ((link->ft == PETSC_VTK_POINT_FIELD) || (link->ft == PETSC_VTK_POINT_VECTOR_FIELD)) hasPoint = PETSC_TRUE;
492:     if ((link->ft == PETSC_VTK_CELL_FIELD) || (link->ft == PETSC_VTK_CELL_VECTOR_FIELD)) hasCell = PETSC_TRUE;
493:   }
494:   if (hasPoint) {
495:     PetscFPrintf(comm, fp, "POINT_DATA %" PetscInt_FMT "\n", totVertices);
496:     for (link = vtk->link; link; link = link->next) {
497:       Vec          X       = (Vec)link->vec;
498:       PetscSection section = NULL, globalSection, newSection = NULL;
499:       char         namebuf[256];
500:       const char  *name;
501:       PetscInt     enforceDof = PETSC_DETERMINE;

503:       if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
504:       if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) enforceDof = 3;
505:       PetscObjectGetName(link->vec, &name);
506:       PetscObjectQuery(link->vec, "section", (PetscObject *)&section);
507:       if (!section) {
508:         DM dmX;

510:         VecGetDM(X, &dmX);
511:         if (dmX) {
512:           DMLabel  subpointMap, subpointMapX;
513:           PetscInt dim, dimX, pStart, pEnd, qStart, qEnd;

515:           DMGetLocalSection(dmX, &section);
516:           /* Here is where we check whether dmX is a submesh of dm */
517:           DMGetDimension(dm, &dim);
518:           DMGetDimension(dmX, &dimX);
519:           DMPlexGetChart(dm, &pStart, &pEnd);
520:           DMPlexGetChart(dmX, &qStart, &qEnd);
521:           DMPlexGetSubpointMap(dm, &subpointMap);
522:           DMPlexGetSubpointMap(dmX, &subpointMapX);
523:           if (((dim != dimX) || ((pEnd - pStart) < (qEnd - qStart))) && subpointMap && !subpointMapX) {
524:             const PetscInt *ind = NULL;
525:             IS              subpointIS;
526:             PetscInt        n = 0, q;

528:             PetscSectionGetChart(section, &qStart, &qEnd);
529:             DMPlexGetSubpointIS(dm, &subpointIS);
530:             if (subpointIS) {
531:               ISGetLocalSize(subpointIS, &n);
532:               ISGetIndices(subpointIS, &ind);
533:             }
534:             PetscSectionCreate(comm, &newSection);
535:             PetscSectionSetChart(newSection, pStart, pEnd);
536:             for (q = qStart; q < qEnd; ++q) {
537:               PetscInt dof, off, p;

539:               PetscSectionGetDof(section, q, &dof);
540:               if (dof) {
541:                 PetscFindInt(q, n, ind, &p);
542:                 if (p >= pStart) {
543:                   PetscSectionSetDof(newSection, p, dof);
544:                   PetscSectionGetOffset(section, q, &off);
545:                   PetscSectionSetOffset(newSection, p, off);
546:                 }
547:               }
548:             }
549:             if (subpointIS) ISRestoreIndices(subpointIS, &ind);
550:             /* No need to setup section */
551:             section = newSection;
552:           }
553:         }
554:       }
556:       if (link->field >= 0) {
557:         const char *fieldname;

559:         PetscSectionGetFieldName(section, link->field, &fieldname);
560:         PetscSectionGetField(section, link->field, &section);
561:         if (fieldname) {
562:           PetscSNPrintf(namebuf, sizeof(namebuf), "%s%s", name, fieldname);
563:         } else {
564:           PetscSNPrintf(namebuf, sizeof(namebuf), "%s%" PetscInt_FMT, name, link->field);
565:         }
566:       } else {
567:         PetscSNPrintf(namebuf, sizeof(namebuf), "%s", name);
568:       }
569:       PetscViewerVTKSanitizeName_Internal(namebuf, sizeof(namebuf));
570:       PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalSection);
571:       for (l = 0; l < loops_per_scalar; l++) DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, namebuf, fp, enforceDof, PETSC_DETERMINE, 1.0, writeComplex, l);
572:       PetscSectionDestroy(&globalSection);
573:       if (newSection) PetscSectionDestroy(&newSection);
574:     }
575:   }
576:   /* Cell Fields */
577:   PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_view_partition", &writePartition, NULL);
578:   if (hasCell || writePartition) {
579:     PetscFPrintf(comm, fp, "CELL_DATA %" PetscInt_FMT "\n", totCells);
580:     for (link = vtk->link; link; link = link->next) {
581:       Vec          X       = (Vec)link->vec;
582:       PetscSection section = NULL, globalSection;
583:       const char  *name    = "";
584:       char         namebuf[256];
585:       PetscInt     enforceDof = PETSC_DETERMINE;

587:       if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
588:       if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) enforceDof = 3;
589:       PetscObjectGetName(link->vec, &name);
590:       PetscObjectQuery(link->vec, "section", (PetscObject *)&section);
591:       if (!section) {
592:         DM dmX;

594:         VecGetDM(X, &dmX);
595:         if (dmX) DMGetLocalSection(dmX, &section);
596:       }
598:       if (link->field >= 0) {
599:         const char *fieldname;

601:         PetscSectionGetFieldName(section, link->field, &fieldname);
602:         PetscSectionGetField(section, link->field, &section);
603:         if (fieldname) {
604:           PetscSNPrintf(namebuf, sizeof(namebuf), "%s%s", name, fieldname);
605:         } else {
606:           PetscSNPrintf(namebuf, sizeof(namebuf), "%s%" PetscInt_FMT, name, link->field);
607:         }
608:       } else {
609:         PetscSNPrintf(namebuf, sizeof(namebuf), "%s", name);
610:       }
611:       PetscViewerVTKSanitizeName_Internal(namebuf, sizeof(namebuf));
612:       PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalSection);
613:       for (l = 0; l < loops_per_scalar; l++) DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, namebuf, fp, enforceDof, PETSC_DETERMINE, 1.0, writeComplex, l);
614:       PetscSectionDestroy(&globalSection);
615:     }
616:     if (writePartition) {
617:       PetscFPrintf(comm, fp, "SCALARS partition int 1\n");
618:       PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n");
619:       DMPlexVTKWritePartition_ASCII(dm, fp);
620:     }
621:   }
622:   /* Cleanup */
623:   PetscFClose(comm, fp);
624:   return 0;
625: }

627: /*@C
628:   DMPlexVTKWriteAll - Write a file containing all the fields that have been provided to the viewer

630:   Collective

632:   Input Parameters:
633: + odm - The `DMPLEX` specifying the mesh, passed as a `PetscObject`
634: - viewer - viewer of type `PETSCVIEWERVTK`

636:   Level: developer

638:   Note:
639:   This function is a callback used by the VTK viewer to actually write the file.
640:   The reason for this odd model is that the VTK file format does not provide any way to write one field at a time.
641:   Instead, metadata for the entire file needs to be available up-front before you can start writing the file.

643: .seealso: [](chapter_unstructured), `DM`, `PETSCVIEWEREXODUSII`, `DMPLEX`, `PETSCVIEWERVTK`
644: @*/
645: PetscErrorCode DMPlexVTKWriteAll(PetscObject odm, PetscViewer viewer)
646: {
647:   DM dm = (DM)odm;

651:   switch (viewer->format) {
652:   case PETSC_VIEWER_ASCII_VTK_DEPRECATED:
653:     DMPlexVTKWriteAll_ASCII(dm, viewer);
654:     break;
655:   case PETSC_VIEWER_VTK_VTU:
656:     DMPlexVTKWriteAll_VTU(dm, viewer);
657:     break;
658:   default:
659:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
660:   }
661:   return 0;
662: }