Actual source code: plexvtk.c

petsc-3.6.1 2015-08-06
Report Typos and Errors
  1: #define PETSCDM_DLL
  2: #include <petsc/private/dmpleximpl.h>    /*I   "petscdmplex.h"   I*/
  3: #include <../src/sys/classes/viewer/impls/vtk/vtkvimpl.h>

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

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

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

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

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

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

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

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

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

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

213:     for (c = 0; c < numCells; ++c) {
214:       DMPlexVTKGetCellType(dm, dim, corners[c], &cellType);
215:       PetscFPrintf(comm, fp, "%d\n", cellType);
216:     }
217:     for (proc = 1; proc < numProcs; ++proc) {
218:       MPI_Status status;

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

238: PetscErrorCode DMPlexVTKWritePartition_ASCII(DM dm, FILE *fp)
239: {
240:   MPI_Comm       comm;
241:   PetscInt       numCells = 0, cellHeight;
242:   PetscInt       numLabelCells, cMax, cStart, cEnd, c;
243:   PetscMPIInt    numProcs, rank, proc, tag;
244:   PetscBool      hasLabel;

248:   PetscObjectGetComm((PetscObject)dm,&comm);
249:   PetscCommGetNewTag(comm, &tag);
250:   MPI_Comm_size(comm, &numProcs);
251:   MPI_Comm_rank(comm, &rank);
252:   DMPlexGetVTKCellHeight(dm, &cellHeight);
253:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
254:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
255:   if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
256:   DMPlexGetStratumSize(dm, "vtk", 1, &numLabelCells);
257:   hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;
258:   for (c = cStart; c < cEnd; ++c) {
259:     if (hasLabel) {
260:       PetscInt value;

262:       DMPlexGetLabelValue(dm, "vtk", c, &value);
263:       if (value != 1) continue;
264:     }
265:     ++numCells;
266:   }
267:   if (!rank) {
268:     for (c = 0; c < numCells; ++c) {PetscFPrintf(comm, fp, "%d\n", rank);}
269:     for (proc = 1; proc < numProcs; ++proc) {
270:       MPI_Status status;

272:       MPI_Recv(&numCells, 1, MPIU_INT, proc, tag, comm, &status);
273:       for (c = 0; c < numCells; ++c) {PetscFPrintf(comm, fp, "%d\n", proc);}
274:     }
275:   } else {
276:     MPI_Send(&numCells, 1, MPIU_INT, 0, tag, comm);
277:   }
278:   return(0);
279: }

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

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

321:       if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
322:           ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
323:         DMPlexGetLabelValue(dm, "vtk", p, &value);
324:         if (value != 1) continue;
325:       }
326:     }
327:     PetscSectionGetDof(section, p, &numDof);
328:     if (numDof) break;
329:   }
330:   MPI_Allreduce(&numDof, &maxDof, 1, MPIU_INT, MPI_MAX, comm);
331:   enforceDof = PetscMax(enforceDof, maxDof);
332:   VecGetArray(v, &array);
333:   if (!rank) {
334:     char formatString[8];

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

629:   Collective

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

635:   Level: developer

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

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

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