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

 73: static PetscErrorCode DMPlexVTKWriteCells_ASCII(DM dm, FILE *fp, PetscInt *totalCells)
 74: {
 75:   MPI_Comm       comm;
 76:   DMLabel        label;
 77:   IS             globalVertexNumbers = NULL;
 78:   const PetscInt *gvertex;
 79:   PetscInt       dim;
 80:   PetscInt       numCorners = 0, totCorners = 0, maxCorners, *corners;
 81:   PetscInt       numCells   = 0, totCells   = 0, maxCells, cellHeight;
 82:   PetscInt       numLabelCells, maxLabelCells, cStart, cEnd, c, vStart, vEnd, v;
 83:   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 %D %D\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, "%D ", nC);
146:       for (v = 0; v < nC; ++v) {
147:         PetscFPrintf(comm, fp, " %D", vertices[v]);
148:       }
149:       PetscFPrintf(comm, fp, "\n");
150:     }
151:     if (size > 1) {PetscMalloc1(maxCorners+maxCells, &remoteVertices);}
152:     for (proc = 1; proc < size; ++proc) {
153:       MPI_Status status;

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

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

175:     PetscMalloc1(numSend, &localVertices);
176:     for (c = cStart, numCells = 0; c < cEnd; ++c) {
177:       PetscInt *closure = NULL;
178:       PetscInt closureSize, value, nC = 0;

180:       if (label) {
181:         DMLabelGetValue(label, c, &value);
182:         if (value != 1) continue;
183:       }
184:       DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
185:       for (v = 0; v < closureSize*2; v += 2) {
186:         if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
187:           const PetscInt gv = gvertex[closure[v] - vStart];
188:           closure[nC++] = gv < 0 ? -(gv+1) : gv;
189:         }
190:       }
191:       corners[numCells++] = nC;
192:       localVertices[k++]  = nC;
193:       for (v = 0; v < nC; ++v, ++k) {
194:         localVertices[k] = closure[v];
195:       }
196:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
197:       DMPlexReorderCell(dm, c, localVertices+k-nC);
198:     }
199:     if (k != numSend) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Invalid number of vertices to send %D should be %D", k, numSend);
200:     MPI_Send(&numSend, 1, MPIU_INT, 0, tag, comm);
201:     MPI_Send(localVertices, numSend, MPIU_INT, 0, tag, comm);
202:     PetscFree(localVertices);
203:   }
204:   ISRestoreIndices(globalVertexNumbers, &gvertex);
205:   PetscFPrintf(comm, fp, "CELL_TYPES %D\n", totCells);
206:   if (rank == 0) {
207:     PetscInt cellType;

209:     for (c = 0; c < numCells; ++c) {
210:       DMPlexVTKGetCellType_Internal(dm, dim, corners[c], &cellType);
211:       PetscFPrintf(comm, fp, "%D\n", cellType);
212:     }
213:     for (proc = 1; proc < size; ++proc) {
214:       MPI_Status status;

216:       MPI_Recv(&numCells, 1, MPIU_INT, proc, tag, comm, &status);
217:       MPI_Recv(corners, numCells, MPIU_INT, proc, tag, comm, &status);
218:       for (c = 0; c < numCells; ++c) {
219:         DMPlexVTKGetCellType_Internal(dm, dim, corners[c], &cellType);
220:         PetscFPrintf(comm, fp, "%D\n", cellType);
221:       }
222:     }
223:   } else {
224:     MPI_Send(&numCells, 1, MPIU_INT, 0, tag, comm);
225:     MPI_Send(corners, numCells, MPIU_INT, 0, tag, comm);
226:   }
227:   PetscFree(corners);
228:   *totalCells = totCells;
229:   return(0);
230: }

232: static PetscErrorCode DMPlexVTKWritePartition_ASCII(DM dm, FILE *fp)
233: {
234:   MPI_Comm       comm;
235:   PetscInt       numCells = 0, cellHeight;
236:   PetscInt       numLabelCells, cStart, cEnd, c;
237:   PetscMPIInt    size, rank, proc, tag;
238:   PetscBool      hasLabel;

242:   PetscObjectGetComm((PetscObject)dm,&comm);
243:   PetscCommGetNewTag(comm, &tag);
244:   MPI_Comm_size(comm, &size);
245:   MPI_Comm_rank(comm, &rank);
246:   DMPlexGetVTKCellHeight(dm, &cellHeight);
247:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
248:   DMGetStratumSize(dm, "vtk", 1, &numLabelCells);
249:   hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;
250:   for (c = cStart; c < cEnd; ++c) {
251:     if (hasLabel) {
252:       PetscInt value;

254:       DMGetLabelValue(dm, "vtk", c, &value);
255:       if (value != 1) continue;
256:     }
257:     ++numCells;
258:   }
259:   if (rank == 0) {
260:     for (c = 0; c < numCells; ++c) {PetscFPrintf(comm, fp, "%d\n", rank);}
261:     for (proc = 1; proc < size; ++proc) {
262:       MPI_Status status;

264:       MPI_Recv(&numCells, 1, MPIU_INT, proc, tag, comm, &status);
265:       for (c = 0; c < numCells; ++c) {PetscFPrintf(comm, fp, "%d\n", proc);}
266:     }
267:   } else {
268:     MPI_Send(&numCells, 1, MPIU_INT, 0, tag, comm);
269:   }
270:   return(0);
271: }

273: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
274: typedef double PetscVTKReal;
275: #elif defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL___FP16)
276: typedef float PetscVTKReal;
277: #else
278: typedef PetscReal PetscVTKReal;
279: #endif

281: static PetscErrorCode DMPlexVTKWriteSection_ASCII(DM dm, PetscSection section, PetscSection globalSection, Vec v, FILE *fp, PetscInt enforceDof, PetscInt precision, PetscReal scale, PetscInt imag)
282: {
283:   MPI_Comm           comm;
284:   const MPI_Datatype mpiType = MPIU_SCALAR;
285:   PetscScalar        *array;
286:   PetscInt           numDof = 0, maxDof;
287:   PetscInt           numLabelCells, cellHeight, cStart, cEnd, numLabelVertices, vStart, vEnd, pStart, pEnd, p;
288:   PetscMPIInt        size, rank, proc, tag;
289:   PetscBool          hasLabel;
290:   PetscErrorCode     ierr;

293:   PetscObjectGetComm((PetscObject)dm,&comm);
296:   if (precision < 0) precision = 6;
297:   PetscCommGetNewTag(comm, &tag);
298:   MPI_Comm_size(comm, &size);
299:   MPI_Comm_rank(comm, &rank);
300:   PetscSectionGetChart(section, &pStart, &pEnd);
301:   /* VTK only wants the values at cells or vertices */
302:   DMPlexGetVTKCellHeight(dm, &cellHeight);
303:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
304:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
305:   pStart   = PetscMax(PetscMin(cStart, vStart), pStart);
306:   pEnd     = PetscMin(PetscMax(cEnd,   vEnd),   pEnd);
307:   DMGetStratumSize(dm, "vtk", 1, &numLabelCells);
308:   DMGetStratumSize(dm, "vtk", 2, &numLabelVertices);
309:   hasLabel = numLabelCells > 0 || numLabelVertices > 0 ? PETSC_TRUE : PETSC_FALSE;
310:   for (p = pStart; p < pEnd; ++p) {
311:     /* Reject points not either cells or vertices */
312:     if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
313:     if (hasLabel) {
314:       PetscInt value;

316:       if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
317:           ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
318:         DMGetLabelValue(dm, "vtk", p, &value);
319:         if (value != 1) continue;
320:       }
321:     }
322:     PetscSectionGetDof(section, p, &numDof);
323:     if (numDof) break;
324:   }
325:   MPIU_Allreduce(&numDof, &maxDof, 1, MPIU_INT, MPI_MAX, comm);
326:   enforceDof = PetscMax(enforceDof, maxDof);
327:   VecGetArray(v, &array);
328:   if (rank == 0) {
329:     PetscVTKReal dval;
330:     PetscScalar  val;
331:     char formatString[8];

333:     PetscSNPrintf(formatString, 8, "%%.%de", precision);
334:     for (p = pStart; p < pEnd; ++p) {
335:       /* Here we lose a way to filter points by keeping them out of the Numbering */
336:       PetscInt dof, off, goff, d;

338:       /* Reject points not either cells or vertices */
339:       if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
340:       if (hasLabel) {
341:         PetscInt value;

343:         if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
344:             ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
345:           DMGetLabelValue(dm, "vtk", p, &value);
346:           if (value != 1) continue;
347:         }
348:       }
349:       PetscSectionGetDof(section, p, &dof);
350:       PetscSectionGetOffset(section, p, &off);
351:       PetscSectionGetOffset(globalSection, p, &goff);
352:       if (dof && goff >= 0) {
353:         for (d = 0; d < dof; d++) {
354:           if (d > 0) {
355:             PetscFPrintf(comm, fp, " ");
356:           }
357:           val = array[off+d];
358:           dval = (PetscVTKReal) ((imag ? PetscImaginaryPart(val) : PetscRealPart(val)) * scale);
359:           PetscFPrintf(comm, fp, formatString, dval);
360:         }
361:         for (d = dof; d < enforceDof; d++) {
362:           PetscFPrintf(comm, fp, " 0.0");
363:         }
364:         PetscFPrintf(comm, fp, "\n");
365:       }
366:     }
367:     for (proc = 1; proc < size; ++proc) {
368:       PetscScalar *remoteValues;
369:       PetscInt    size = 0, d;
370:       MPI_Status  status;

372:       MPI_Recv(&size, 1, MPIU_INT, proc, tag, comm, &status);
373:       PetscMalloc1(size, &remoteValues);
374:       MPI_Recv(remoteValues, size, mpiType, proc, tag, comm, &status);
375:       for (p = 0; p < size/maxDof; ++p) {
376:         for (d = 0; d < maxDof; ++d) {
377:           if (d > 0) {
378:             PetscFPrintf(comm, fp, " ");
379:           }
380:           val = remoteValues[p*maxDof+d];
381:           dval = (PetscVTKReal) ((imag ? PetscImaginaryPart(val) : PetscRealPart(val)) * scale);
382:           PetscFPrintf(comm, fp, formatString, dval);
383:         }
384:         for (d = maxDof; d < enforceDof; ++d) {
385:           PetscFPrintf(comm, fp, " 0.0");
386:         }
387:         PetscFPrintf(comm, fp, "\n");
388:       }
389:       PetscFree(remoteValues);
390:     }
391:   } else {
392:     PetscScalar *localValues;
393:     PetscInt    size, k = 0;

395:     PetscSectionGetStorageSize(section, &size);
396:     PetscMalloc1(size, &localValues);
397:     for (p = pStart; p < pEnd; ++p) {
398:       PetscInt dof, off, goff, d;

400:       /* Reject points not either cells or vertices */
401:       if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
402:       if (hasLabel) {
403:         PetscInt value;

405:         if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
406:             ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
407:           DMGetLabelValue(dm, "vtk", p, &value);
408:           if (value != 1) continue;
409:         }
410:       }
411:       PetscSectionGetDof(section, p, &dof);
412:       PetscSectionGetOffset(section, p, &off);
413:       PetscSectionGetOffset(globalSection, p, &goff);
414:       if (goff >= 0) {
415:         for (d = 0; d < dof; ++d) {
416:           localValues[k++] = array[off+d];
417:         }
418:       }
419:     }
420:     MPI_Send(&k, 1, MPIU_INT, 0, tag, comm);
421:     MPI_Send(localValues, k, mpiType, 0, tag, comm);
422:     PetscFree(localValues);
423:   }
424:   VecRestoreArray(v, &array);
425:   return(0);
426: }

428: 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)
429: {
430:   MPI_Comm       comm;
431:   PetscInt       numDof = 0, maxDof;
432:   PetscInt       pStart, pEnd, p;

436:   PetscObjectGetComm((PetscObject)dm,&comm);
437:   PetscSectionGetChart(section, &pStart, &pEnd);
438:   for (p = pStart; p < pEnd; ++p) {
439:     PetscSectionGetDof(section, p, &numDof);
440:     if (numDof) break;
441:   }
442:   numDof = PetscMax(numDof, enforceDof);
443:   MPIU_Allreduce(&numDof, &maxDof, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
444:   if (!name) name = "Unknown";
445:   if (maxDof == 3) {
446:     if (nameComplex) {
447:       PetscFPrintf(comm, fp, "VECTORS %s.%s double\n", name, imag ? "Im" : "Re");
448:     } else {
449:       PetscFPrintf(comm, fp, "VECTORS %s double\n", name);
450:     }
451:   } else {
452:     if (nameComplex) {
453:       PetscFPrintf(comm, fp, "SCALARS %s.%s double %D\n", name, imag ? "Im" : "Re", maxDof);
454:     } else {
455:       PetscFPrintf(comm, fp, "SCALARS %s double %D\n", name, maxDof);
456:     }
457:     PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n");
458:   }
459:   DMPlexVTKWriteSection_ASCII(dm, section, globalSection, field, fp, enforceDof, precision, scale, imag);
460:   return(0);
461: }

463: static PetscErrorCode DMPlexVTKWriteAll_ASCII(DM dm, PetscViewer viewer)
464: {
465:   MPI_Comm                 comm;
466:   PetscViewer_VTK          *vtk = (PetscViewer_VTK*) viewer->data;
467:   FILE                     *fp;
468:   PetscViewerVTKObjectLink link;
469:   PetscSection             coordSection, globalCoordSection;
470:   PetscLayout              vLayout;
471:   Vec                      coordinates;
472:   PetscReal                lengthScale;
473:   PetscInt                 totVertices, totCells = 0, loops_per_scalar, l;
474:   PetscBool                hasPoint = PETSC_FALSE, hasCell = PETSC_FALSE, writePartition = PETSC_FALSE, localized, writeComplex;
475:   const char               *dmname;
476:   PetscErrorCode           ierr;

479: #if defined(PETSC_USE_COMPLEX)
480:   loops_per_scalar = 2;
481:   writeComplex = PETSC_TRUE;
482: #else
483:   loops_per_scalar = 1;
484:   writeComplex = PETSC_FALSE;
485: #endif
486:   DMGetCoordinatesLocalized(dm,&localized);
487:   if (localized) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"VTK output with localized coordinates not yet supported");
488:   PetscObjectGetComm((PetscObject)dm,&comm);
489:   PetscFOpen(comm, vtk->filename, "wb", &fp);
490:   PetscObjectGetName((PetscObject)dm, &dmname);
491:   PetscFPrintf(comm, fp, "# vtk DataFile Version 2.0\n");
492:   PetscFPrintf(comm, fp, "%s\n", dmname);
493:   PetscFPrintf(comm, fp, "ASCII\n");
494:   PetscFPrintf(comm, fp, "DATASET UNSTRUCTURED_GRID\n");
495:   /* Vertices */
496:   DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
497:   DMGetCoordinateSection(dm, &coordSection);
498:   PetscSectionCreateGlobalSection(coordSection, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalCoordSection);
499:   DMGetCoordinatesLocal(dm, &coordinates);
500:   PetscSectionGetPointLayout(PetscObjectComm((PetscObject)dm), globalCoordSection, &vLayout);
501:   PetscLayoutGetSize(vLayout, &totVertices);
502:   PetscFPrintf(comm, fp, "POINTS %D double\n", totVertices);
503:   DMPlexVTKWriteSection_ASCII(dm, coordSection, globalCoordSection, coordinates, fp, 3, PETSC_DETERMINE, lengthScale, 0);
504:   /* Cells */
505:   DMPlexVTKWriteCells_ASCII(dm, fp, &totCells);
506:   /* Vertex fields */
507:   for (link = vtk->link; link; link = link->next) {
508:     if ((link->ft == PETSC_VTK_POINT_FIELD) || (link->ft == PETSC_VTK_POINT_VECTOR_FIELD)) hasPoint = PETSC_TRUE;
509:     if ((link->ft == PETSC_VTK_CELL_FIELD)  || (link->ft == PETSC_VTK_CELL_VECTOR_FIELD))  hasCell  = PETSC_TRUE;
510:   }
511:   if (hasPoint) {
512:     PetscFPrintf(comm, fp, "POINT_DATA %D\n", totVertices);
513:     for (link = vtk->link; link; link = link->next) {
514:       Vec          X = (Vec) link->vec;
515:       PetscSection section = NULL, globalSection, newSection = NULL;
516:       char         namebuf[256];
517:       const char   *name;
518:       PetscInt     enforceDof = PETSC_DETERMINE;

520:       if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
521:       if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) enforceDof = 3;
522:       PetscObjectGetName(link->vec, &name);
523:       PetscObjectQuery(link->vec, "section", (PetscObject*) &section);
524:       if (!section) {
525:         DM           dmX;

527:         VecGetDM(X, &dmX);
528:         if (dmX) {
529:           DMLabel  subpointMap, subpointMapX;
530:           PetscInt dim, dimX, pStart, pEnd, qStart, qEnd;

532:           DMGetLocalSection(dmX, &section);
533:           /* Here is where we check whether dmX is a submesh of dm */
534:           DMGetDimension(dm,  &dim);
535:           DMGetDimension(dmX, &dimX);
536:           DMPlexGetChart(dm,  &pStart, &pEnd);
537:           DMPlexGetChart(dmX, &qStart, &qEnd);
538:           DMPlexGetSubpointMap(dm,  &subpointMap);
539:           DMPlexGetSubpointMap(dmX, &subpointMapX);
540:           if (((dim != dimX) || ((pEnd-pStart) < (qEnd-qStart))) && subpointMap && !subpointMapX) {
541:             const PetscInt *ind = NULL;
542:             IS              subpointIS;
543:             PetscInt        n = 0, q;

545:             PetscSectionGetChart(section, &qStart, &qEnd);
546:             DMPlexGetSubpointIS(dm, &subpointIS);
547:             if (subpointIS) {
548:               ISGetLocalSize(subpointIS, &n);
549:               ISGetIndices(subpointIS, &ind);
550:             }
551:             PetscSectionCreate(comm, &newSection);
552:             PetscSectionSetChart(newSection, pStart, pEnd);
553:             for (q = qStart; q < qEnd; ++q) {
554:               PetscInt dof, off, p;

556:               PetscSectionGetDof(section, q, &dof);
557:               if (dof) {
558:                 PetscFindInt(q, n, ind, &p);
559:                 if (p >= pStart) {
560:                   PetscSectionSetDof(newSection, p, dof);
561:                   PetscSectionGetOffset(section, q, &off);
562:                   PetscSectionSetOffset(newSection, p, off);
563:                 }
564:               }
565:             }
566:             if (subpointIS) {
567:               ISRestoreIndices(subpointIS, &ind);
568:             }
569:             /* No need to setup section */
570:             section = newSection;
571:           }
572:         }
573:       }
574:       if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it and could not create one from VecGetDM()", name);
575:       if (link->field >= 0) {
576:         const char *fieldname;

578:         PetscSectionGetFieldName(section, link->field, &fieldname);
579:         PetscSectionGetField(section, link->field, &section);
580:         if (fieldname) {
581:           PetscSNPrintf(namebuf, sizeof(namebuf), "%s%s", name, fieldname);
582:         } else {
583:           PetscSNPrintf(namebuf, sizeof(namebuf), "%s%D", name, link->field);
584:         }
585:       } else {
586:         PetscSNPrintf(namebuf, sizeof(namebuf), "%s", name);
587:       }
588:       PetscViewerVTKSanitizeName_Internal(namebuf, sizeof(namebuf));
589:       PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalSection);
590:       for (l = 0; l < loops_per_scalar; l++) {
591:         DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, namebuf, fp, enforceDof, PETSC_DETERMINE, 1.0, writeComplex, l);
592:       }
593:       PetscSectionDestroy(&globalSection);
594:       if (newSection) {PetscSectionDestroy(&newSection);}
595:     }
596:   }
597:   /* Cell Fields */
598:   PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_view_partition", &writePartition, NULL);
599:   if (hasCell || writePartition) {
600:     PetscFPrintf(comm, fp, "CELL_DATA %D\n", totCells);
601:     for (link = vtk->link; link; link = link->next) {
602:       Vec          X = (Vec) link->vec;
603:       PetscSection section = NULL, globalSection;
604:       const char   *name = "";
605:       char         namebuf[256];
606:       PetscInt     enforceDof = PETSC_DETERMINE;

608:       if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
609:       if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) enforceDof = 3;
610:       PetscObjectGetName(link->vec, &name);
611:       PetscObjectQuery(link->vec, "section", (PetscObject*) &section);
612:       if (!section) {
613:         DM           dmX;

615:         VecGetDM(X, &dmX);
616:         if (dmX) {
617:           DMGetLocalSection(dmX, &section);
618:         }
619:       }
620:       if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it and could not create one from VecGetDM()", name);
621:       if (link->field >= 0) {
622:         const char *fieldname;

624:         PetscSectionGetFieldName(section, link->field, &fieldname);
625:         PetscSectionGetField(section, link->field, &section);
626:         if (fieldname) {
627:           PetscSNPrintf(namebuf, sizeof(namebuf), "%s%s", name, fieldname);
628:         } else {
629:           PetscSNPrintf(namebuf, sizeof(namebuf), "%s%D", name, link->field);
630:         }
631:       } else {
632:         PetscSNPrintf(namebuf, sizeof(namebuf), "%s", name);
633:       }
634:       PetscViewerVTKSanitizeName_Internal(namebuf, sizeof(namebuf));
635:       PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalSection);
636:       for (l = 0; l < loops_per_scalar; l++) {
637:         DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, namebuf, fp, enforceDof, PETSC_DETERMINE, 1.0, writeComplex, l);
638:       }
639:       PetscSectionDestroy(&globalSection);
640:     }
641:     if (writePartition) {
642:       PetscFPrintf(comm, fp, "SCALARS partition int 1\n");
643:       PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n");
644:       DMPlexVTKWritePartition_ASCII(dm, fp);
645:     }
646:   }
647:   /* Cleanup */
648:   PetscSectionDestroy(&globalCoordSection);
649:   PetscLayoutDestroy(&vLayout);
650:   PetscFClose(comm, fp);
651:   return(0);
652: }

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

657:   Collective

659:   Input Parameters:
660: + odm - The DMPlex specifying the mesh, passed as a PetscObject
661: - viewer - viewer of type VTK

663:   Level: developer

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

670: .seealso: PETSCVIEWERVTK
671: @*/
672: PetscErrorCode DMPlexVTKWriteAll(PetscObject odm, PetscViewer viewer)
673: {
674:   DM             dm = (DM) odm;
675:   PetscBool      isvtk;

681:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);
682:   if (!isvtk) SETERRQ1(PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_INCOMP, "Cannot use viewer type %s", ((PetscObject)viewer)->type_name);
683:   switch (viewer->format) {
684:   case PETSC_VIEWER_ASCII_VTK_DEPRECATED:
685:     DMPlexVTKWriteAll_ASCII(dm, viewer);
686:     break;
687:   case PETSC_VIEWER_VTK_VTU:
688:     DMPlexVTKWriteAll_VTU(dm, viewer);
689:     break;
690:   default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
691:   }
692:   return(0);
693: }