Actual source code: plexvtk.c

petsc-3.9.4 2018-09-11
Report Typos and Errors
  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 8:
 55:       *cellType = 12; /* VTK_HEXAHEDRON */
 56:       break;
 57:     case 10:
 58:       *cellType = 24; /* VTK_QUADRATIC_TETRA */
 59:       break;
 60:     case 27:
 61:       *cellType = 29; /* VTK_QUADRATIC_HEXAHEDRON */
 62:       break;
 63:     default:
 64:       break;
 65:     }
 66:   }
 67:   return(0);
 68: }

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

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

102:     if (label) {
103:       DMLabelGetValue(label, c, &value);
104:       if (value != 1) continue;
105:     }
106:     DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
107:     for (v = 0; v < closureSize*2; v += 2) {
108:       if ((closure[v] >= vStart) && (closure[v] < vEnd)) ++numCorners;
109:     }
110:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
111:     ++numCells;
112:   }
113:   maxCells = numCells;
114:   MPI_Reduce(&numCells, &totCells, 1, MPIU_INT, MPI_SUM, 0, comm);
115:   MPI_Reduce(&numCells, &maxCells, 1, MPIU_INT, MPI_MAX, 0, comm);
116:   MPI_Reduce(&numCorners, &totCorners, 1, MPIU_INT, MPI_SUM, 0, comm);
117:   MPI_Reduce(&numCorners, &maxCorners, 1, MPIU_INT, MPI_MAX, 0, comm);
118:   DMPlexGetVertexNumbering(dm, &globalVertexNumbers);
119:   ISGetIndices(globalVertexNumbers, &gvertex);
120:   PetscMalloc1(maxCells, &corners);
121:   PetscFPrintf(comm, fp, "CELLS %D %D\n", totCells, totCorners+totCells);
122:   if (!rank) {
123:     PetscInt *remoteVertices;
124:     int      *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:       corners[numCells++] = nC;
144:       PetscFPrintf(comm, fp, "%D ", nC);
145:       DMPlexInvertCell(dim, nC, vertices);
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:         DMPlexInvertCell(dim, nC, vertices);
164:         PetscFPrintf(comm, fp, "%D ", nC);
165:         for (v = 0; v < nC; ++v) {
166:           PetscFPrintf(comm, fp, " %D", vertices[v]);
167:         }
168:         PetscFPrintf(comm, fp, "\n");
169:       }
170:     }
171:     if (size > 1) {PetscFree(remoteVertices);}
172:     PetscFree(vertices);
173:   } else {
174:     PetscInt *localVertices, numSend = numCells+numCorners, k = 0;

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

181:       if (label) {
182:         DMLabelGetValue(label, c, &value);
183:         if (value != 1) continue;
184:       }
185:       DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
186:       for (v = 0; v < closureSize*2; v += 2) {
187:         if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
188:           const PetscInt gv = gvertex[closure[v] - vStart];
189:           closure[nC++] = gv < 0 ? -(gv+1) : gv;
190:         }
191:       }
192:       corners[numCells++] = nC;
193:       localVertices[k++]  = nC;
194:       for (v = 0; v < nC; ++v, ++k) {
195:         localVertices[k] = closure[v];
196:       }
197:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
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) {
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, cMax, 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:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
249:   if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
250:   DMGetStratumSize(dm, "vtk", 1, &numLabelCells);
251:   hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;
252:   for (c = cStart; c < cEnd; ++c) {
253:     if (hasLabel) {
254:       PetscInt value;

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

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

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

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

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

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

333:       /* Reject points not either cells or vertices */
334:       if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
335:       if (hasLabel) {
336:         PetscInt value;

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

365:       MPI_Recv(&size, 1, MPIU_INT, proc, tag, comm, &status);
366:       PetscMalloc1(size, &remoteValues);
367:       MPI_Recv(remoteValues, size, mpiType, proc, tag, comm, &status);
368:       for (p = 0; p < size/maxDof; ++p) {
369:         for (d = 0; d < maxDof; ++d) {
370:           if (d > 0) {
371:             PetscFPrintf(comm, fp, " ");
372:           }
373:           PetscFPrintf(comm, fp, formatString, PetscRealPart(remoteValues[p*maxDof+d])*scale);
374:         }
375:         for (d = maxDof; d < enforceDof; ++d) {
376:           PetscFPrintf(comm, fp, " 0.0");
377:         }
378:         PetscFPrintf(comm, fp, "\n");
379:       }
380:       PetscFree(remoteValues);
381:     }
382:   } else {
383:     PetscScalar *localValues;
384:     PetscInt    size, k = 0;

386:     PetscSectionGetStorageSize(section, &size);
387:     PetscMalloc1(size, &localValues);
388:     for (p = pStart; p < pEnd; ++p) {
389:       PetscInt dof, off, goff, d;

391:       /* Reject points not either cells or vertices */
392:       if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
393:       if (hasLabel) {
394:         PetscInt value;

396:         if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
397:             ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
398:           DMGetLabelValue(dm, "vtk", p, &value);
399:           if (value != 1) continue;
400:         }
401:       }
402:       PetscSectionGetDof(section, p, &dof);
403:       PetscSectionGetOffset(section, p, &off);
404:       PetscSectionGetOffset(globalSection, p, &goff);
405:       if (goff >= 0) {
406:         for (d = 0; d < dof; ++d) {
407:           localValues[k++] = array[off+d];
408:         }
409:       }
410:     }
411:     MPI_Send(&k, 1, MPIU_INT, 0, tag, comm);
412:     MPI_Send(localValues, k, mpiType, 0, tag, comm);
413:     PetscFree(localValues);
414:   }
415:   VecRestoreArray(v, &array);
416:   return(0);
417: }

419: static PetscErrorCode DMPlexVTKWriteField_ASCII(DM dm, PetscSection section, PetscSection globalSection, Vec field, const char name[], FILE *fp, PetscInt enforceDof, PetscInt precision, PetscReal scale)
420: {
421:   MPI_Comm       comm;
422:   PetscInt       numDof = 0, maxDof;
423:   PetscInt       pStart, pEnd, p;

427:   PetscObjectGetComm((PetscObject)dm,&comm);
428:   PetscSectionGetChart(section, &pStart, &pEnd);
429:   for (p = pStart; p < pEnd; ++p) {
430:     PetscSectionGetDof(section, p, &numDof);
431:     if (numDof) break;
432:   }
433:   numDof = PetscMax(numDof, enforceDof);
434:   MPIU_Allreduce(&numDof, &maxDof, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
435:   if (!name) name = "Unknown";
436:   if (maxDof == 3) {
437:     PetscFPrintf(comm, fp, "VECTORS %s double\n", name);
438:   } else {
439:     PetscFPrintf(comm, fp, "SCALARS %s double %D\n", name, maxDof);
440:     PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n");
441:   }
442:   DMPlexVTKWriteSection_ASCII(dm, section, globalSection, field, fp, enforceDof, precision, scale);
443:   return(0);
444: }

446: static PetscErrorCode DMPlexVTKWriteAll_ASCII(DM dm, PetscViewer viewer)
447: {
448:   MPI_Comm                 comm;
449:   PetscViewer_VTK          *vtk = (PetscViewer_VTK*) viewer->data;
450:   FILE                     *fp;
451:   PetscViewerVTKObjectLink link;
452:   PetscSection             coordSection, globalCoordSection;
453:   PetscLayout              vLayout;
454:   Vec                      coordinates;
455:   PetscReal                lengthScale;
456:   PetscInt                 vMax, totVertices, totCells = 0;
457:   PetscBool                hasPoint = PETSC_FALSE, hasCell = PETSC_FALSE, writePartition = PETSC_FALSE, localized;
458:   PetscErrorCode           ierr;

461:   DMGetCoordinatesLocalized(dm,&localized);
462:   if (localized) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"VTK output with localized coordinates not yet supported");
463:   PetscObjectGetComm((PetscObject)dm,&comm);
464:   PetscFOpen(comm, vtk->filename, "wb", &fp);
465:   PetscFPrintf(comm, fp, "# vtk DataFile Version 2.0\n");
466:   PetscFPrintf(comm, fp, "Simplicial Mesh Example\n");
467:   PetscFPrintf(comm, fp, "ASCII\n");
468:   PetscFPrintf(comm, fp, "DATASET UNSTRUCTURED_GRID\n");
469:   /* Vertices */
470:   DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
471:   DMGetCoordinateSection(dm, &coordSection);
472:   PetscSectionCreateGlobalSection(coordSection, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalCoordSection);
473:   DMGetCoordinatesLocal(dm, &coordinates);
474:   DMPlexGetHybridBounds(dm, NULL, NULL, NULL, &vMax);
475:   if (vMax >= 0) {
476:     PetscInt pStart, pEnd, p, localSize = 0;

478:     PetscSectionGetChart(globalCoordSection, &pStart, &pEnd);
479:     pEnd = PetscMin(pEnd, vMax);
480:     for (p = pStart; p < pEnd; ++p) {
481:       PetscInt dof;

483:       PetscSectionGetDof(globalCoordSection, p, &dof);
484:       if (dof > 0) ++localSize;
485:     }
486:     PetscLayoutCreate(PetscObjectComm((PetscObject)dm), &vLayout);
487:     PetscLayoutSetLocalSize(vLayout, localSize);
488:     PetscLayoutSetBlockSize(vLayout, 1);
489:     PetscLayoutSetUp(vLayout);
490:   } else {
491:     PetscSectionGetPointLayout(PetscObjectComm((PetscObject)dm), globalCoordSection, &vLayout);
492:   }
493:   PetscLayoutGetSize(vLayout, &totVertices);
494:   PetscFPrintf(comm, fp, "POINTS %D double\n", totVertices);
495:   DMPlexVTKWriteSection_ASCII(dm, coordSection, globalCoordSection, coordinates, fp, 3, PETSC_DETERMINE, lengthScale);
496:   /* Cells */
497:   DMPlexVTKWriteCells_ASCII(dm, fp, &totCells);
498:   /* Vertex fields */
499:   for (link = vtk->link; link; link = link->next) {
500:     if ((link->ft == PETSC_VTK_POINT_FIELD) || (link->ft == PETSC_VTK_POINT_VECTOR_FIELD)) hasPoint = PETSC_TRUE;
501:     if ((link->ft == PETSC_VTK_CELL_FIELD)  || (link->ft == PETSC_VTK_CELL_VECTOR_FIELD))  hasCell  = PETSC_TRUE;
502:   }
503:   if (hasPoint) {
504:     PetscFPrintf(comm, fp, "POINT_DATA %D\n", totVertices);
505:     for (link = vtk->link; link; link = link->next) {
506:       Vec          X = (Vec) link->vec;
507:       DM           dmX;
508:       PetscSection section, globalSection, newSection = NULL;
509:       const char   *name;
510:       PetscInt     enforceDof = PETSC_DETERMINE;

512:       if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
513:       if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) enforceDof = 3;
514:       PetscObjectGetName(link->vec, &name);
515:       VecGetDM(X, &dmX);
516:       if (dmX) {
517:         DMLabel  subpointMap, subpointMapX;
518:         PetscInt dim, dimX, pStart, pEnd, qStart, qEnd;

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

533:           PetscSectionGetChart(section, &qStart, &qEnd);
534:           DMPlexCreateSubpointIS(dm, &subpointIS);
535:           if (subpointIS) {
536:             ISGetLocalSize(subpointIS, &n);
537:             ISGetIndices(subpointIS, &ind);
538:           }
539:           PetscSectionCreate(comm, &newSection);
540:           PetscSectionSetChart(newSection, pStart, pEnd);
541:           for (q = qStart; q < qEnd; ++q) {
542:             PetscInt dof, off, p;

544:             PetscSectionGetDof(section, q, &dof);
545:             if (dof) {
546:               PetscFindInt(q, n, ind, &p);
547:               if (p >= pStart) {
548:                 PetscSectionSetDof(newSection, p, dof);
549:                 PetscSectionGetOffset(section, q, &off);
550:                 PetscSectionSetOffset(newSection, p, off);
551:               }
552:             }
553:           }
554:           if (subpointIS) {
555:             ISRestoreIndices(subpointIS, &ind);
556:             ISDestroy(&subpointIS);
557:           }
558:           /* No need to setup section */
559:           section = newSection;
560:         }
561:       } else {
562:         PetscObjectQuery(link->vec, "section", (PetscObject*) &section);
563:         if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
564:       }
565:       if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
566:       PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalSection);
567:       DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, name, fp, enforceDof, PETSC_DETERMINE, 1.0);
568:       PetscSectionDestroy(&globalSection);
569:       if (newSection) {PetscSectionDestroy(&newSection);}
570:     }
571:   }
572:   /* Cell Fields */
573:   PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_view_partition", &writePartition, NULL);
574:   if (hasCell || writePartition) {
575:     PetscFPrintf(comm, fp, "CELL_DATA %D\n", totCells);
576:     for (link = vtk->link; link; link = link->next) {
577:       Vec          X = (Vec) link->vec;
578:       DM           dmX;
579:       PetscSection section, globalSection;
580:       const char   *name;
581:       PetscInt     enforceDof = PETSC_DETERMINE;

583:       if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
584:       if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) enforceDof = 3;
585:       PetscObjectGetName(link->vec, &name);
586:       VecGetDM(X, &dmX);
587:       if (dmX) {
588:         DMGetDefaultSection(dmX, &section);
589:       } else {
590:         PetscContainer c;

592:         PetscObjectQuery(link->vec, "section", (PetscObject*) &c);
593:         if (!c) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
594:         PetscContainerGetPointer(c, (void**) &section);
595:       }
596:       if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
597:       PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalSection);
598:       DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, name, fp, enforceDof, PETSC_DETERMINE, 1.0);
599:       PetscSectionDestroy(&globalSection);
600:     }
601:     if (writePartition) {
602:       PetscFPrintf(comm, fp, "SCALARS partition int 1\n");
603:       PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n");
604:       DMPlexVTKWritePartition_ASCII(dm, fp);
605:     }
606:   }
607:   /* Cleanup */
608:   PetscSectionDestroy(&globalCoordSection);
609:   PetscLayoutDestroy(&vLayout);
610:   PetscFClose(comm, fp);
611:   return(0);
612: }

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

617:   Collective

619:   Input Arguments:
620: + odm - The DMPlex specifying the mesh, passed as a PetscObject
621: - viewer - viewer of type VTK

623:   Level: developer

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

630: .seealso: PETSCVIEWERVTK
631: @*/
632: PetscErrorCode DMPlexVTKWriteAll(PetscObject odm, PetscViewer viewer)
633: {
634:   DM             dm = (DM) odm;
635:   PetscBool      isvtk;

641:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);
642:   if (!isvtk) SETERRQ1(PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_INCOMP, "Cannot use viewer type %s", ((PetscObject)viewer)->type_name);
643:   switch (viewer->format) {
644:   case PETSC_VIEWER_ASCII_VTK:
645:     DMPlexVTKWriteAll_ASCII(dm, viewer);
646:     break;
647:   case PETSC_VIEWER_VTK_VTU:
648:     DMPlexVTKWriteAll_VTU(dm, viewer);
649:     break;
650:   default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
651:   }
652:   return(0);
653: }