Actual source code: plexvtk.c

petsc-3.11.4 2019-09-28
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 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) {
124:     PetscInt *remoteVertices;
125:     int      *vertices;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

430: static PetscErrorCode DMPlexVTKWriteField_ASCII(DM dm, PetscSection section, PetscSection globalSection, Vec field, const char name[], FILE *fp, PetscInt enforceDof, PetscInt precision, PetscReal scale)
431: {
432:   MPI_Comm       comm;
433:   PetscInt       numDof = 0, maxDof;
434:   PetscInt       pStart, pEnd, p;

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

457: static PetscErrorCode DMPlexVTKWriteAll_ASCII(DM dm, PetscViewer viewer)
458: {
459:   MPI_Comm                 comm;
460:   PetscViewer_VTK          *vtk = (PetscViewer_VTK*) viewer->data;
461:   FILE                     *fp;
462:   PetscViewerVTKObjectLink link;
463:   PetscSection             coordSection, globalCoordSection;
464:   PetscLayout              vLayout;
465:   Vec                      coordinates;
466:   PetscReal                lengthScale;
467:   PetscInt                 vMax, totVertices, totCells = 0;
468:   PetscBool                hasPoint = PETSC_FALSE, hasCell = PETSC_FALSE, writePartition = PETSC_FALSE, localized;
469:   PetscErrorCode           ierr;

472:   DMGetCoordinatesLocalized(dm,&localized);
473:   if (localized) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"VTK output with localized coordinates not yet supported");
474:   PetscObjectGetComm((PetscObject)dm,&comm);
475:   PetscFOpen(comm, vtk->filename, "wb", &fp);
476:   PetscFPrintf(comm, fp, "# vtk DataFile Version 2.0\n");
477:   PetscFPrintf(comm, fp, "Simplicial Mesh Example\n");
478:   PetscFPrintf(comm, fp, "ASCII\n");
479:   PetscFPrintf(comm, fp, "DATASET UNSTRUCTURED_GRID\n");
480:   /* Vertices */
481:   DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
482:   DMGetCoordinateSection(dm, &coordSection);
483:   PetscSectionCreateGlobalSection(coordSection, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalCoordSection);
484:   DMGetCoordinatesLocal(dm, &coordinates);
485:   DMPlexGetHybridBounds(dm, NULL, NULL, NULL, &vMax);
486:   if (vMax >= 0) {
487:     PetscInt pStart, pEnd, p, localSize = 0;

489:     PetscSectionGetChart(globalCoordSection, &pStart, &pEnd);
490:     pEnd = PetscMin(pEnd, vMax);
491:     for (p = pStart; p < pEnd; ++p) {
492:       PetscInt dof;

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

523:       if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
524:       if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) enforceDof = 3;
525:       PetscObjectGetName(link->vec, &name);
526:       VecGetDM(X, &dmX);
527:       if (dmX) {
528:         DMLabel  subpointMap, subpointMapX;
529:         PetscInt dim, dimX, pStart, pEnd, qStart, qEnd;

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

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

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

594:       if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
595:       if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) enforceDof = 3;
596:       PetscObjectGetName(link->vec, &name);
597:       VecGetDM(X, &dmX);
598:       if (dmX) {
599:         DMGetSection(dmX, &section);
600:       } else {
601:         PetscContainer c;

603:         PetscObjectQuery(link->vec, "section", (PetscObject*) &c);
604:         if (!c) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
605:         PetscContainerGetPointer(c, (void**) &section);
606:       }
607:       if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
608:       PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalSection);
609:       DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, name, fp, enforceDof, PETSC_DETERMINE, 1.0);
610:       PetscSectionDestroy(&globalSection);
611:     }
612:     if (writePartition) {
613:       PetscFPrintf(comm, fp, "SCALARS partition int 1\n");
614:       PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n");
615:       DMPlexVTKWritePartition_ASCII(dm, fp);
616:     }
617:   }
618:   /* Cleanup */
619:   PetscSectionDestroy(&globalCoordSection);
620:   PetscLayoutDestroy(&vLayout);
621:   PetscFClose(comm, fp);
622:   return(0);
623: }

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

628:   Collective

630:   Input Arguments:
631: + odm - The DMPlex specifying the mesh, passed as a PetscObject
632: - viewer - viewer of type VTK

634:   Level: developer

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

641: .seealso: PETSCVIEWERVTK
642: @*/
643: PetscErrorCode DMPlexVTKWriteAll(PetscObject odm, PetscViewer viewer)
644: {
645:   DM             dm = (DM) odm;
646:   PetscBool      isvtk;

652:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);
653:   if (!isvtk) SETERRQ1(PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_INCOMP, "Cannot use viewer type %s", ((PetscObject)viewer)->type_name);
654:   switch (viewer->format) {
655:   case PETSC_VIEWER_ASCII_VTK:
656:     DMPlexVTKWriteAll_ASCII(dm, viewer);
657:     break;
658:   case PETSC_VIEWER_VTK_VTU:
659:     DMPlexVTKWriteAll_VTU(dm, viewer);
660:     break;
661:   default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
662:   }
663:   return(0);
664: }