Actual source code: plexvtk.c

petsc-3.4.4 2014-03-13
  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:   IS             globalVertexNumbers;
 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, cMax, cStart, cEnd, c, vStart, vEnd, v;
 83:   PetscMPIInt    numProcs, rank, proc, tag;
 84:   PetscBool      hasLabel;

 88:   PetscObjectGetComm((PetscObject)dm,&comm);
 89:   PetscCommGetNewTag(comm, &tag);
 90:   MPI_Comm_size(comm, &numProcs);
 91:   MPI_Comm_rank(comm, &rank);
 92:   DMPlexGetDimension(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:   DMPlexGetStratumSize(dm, "vtk", 1, &numLabelCells);

100:   hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;
101:   for (c = cStart; c < cEnd; ++c) {
102:     PetscInt *closure = NULL;
103:     PetscInt closureSize;

105:     if (hasLabel) {
106:       PetscInt value;

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

130:     for (c = cStart, numCells = 0; c < cEnd; ++c) {
131:       PetscInt *closure = NULL;
132:       PetscInt closureSize, nC = 0, tmp;

134:       if (hasLabel) {
135:         PetscInt value;

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

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

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

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

183:       if (hasLabel) {
184:         PetscInt value;

186:         DMPlexGetLabelValue(dm, "vtk", 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 DMPlexVTKWriteSection_ASCII(DM dm, PetscSection section, PetscSection globalSection, Vec v, FILE *fp, PetscInt enforceDof, PetscInt precision, PetscReal scale)
239: {
240:   MPI_Comm           comm;
241:   const MPI_Datatype mpiType = MPIU_SCALAR;
242:   PetscScalar        *array;
243:   PetscInt           numDof = 0, maxDof;
244:   PetscInt           numLabelCells, cellHeight, cMax, cStart, cEnd, numLabelVertices, vMax, vStart, vEnd, pStart, pEnd, p;
245:   PetscMPIInt        numProcs, rank, proc, tag;
246:   PetscBool          hasLabel;
247:   PetscErrorCode     ierr;

250:   PetscObjectGetComm((PetscObject)dm,&comm);
253:   if (precision < 0) precision = 6;
254:   PetscCommGetNewTag(comm, &tag);
255:   MPI_Comm_size(comm, &numProcs);
256:   MPI_Comm_rank(comm, &rank);
257:   PetscSectionGetChart(section, &pStart, &pEnd);
258:   /* VTK only wants the values at cells or vertices */
259:   DMPlexGetVTKCellHeight(dm, &cellHeight);
260:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
261:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
262:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, &vMax);
263:   if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
264:   if (vMax >= 0) vEnd = PetscMin(vEnd, vMax);
265:   pStart   = PetscMax(PetscMin(cStart, vStart), pStart);
266:   pEnd     = PetscMin(PetscMax(cEnd,   vEnd),   pEnd);
267:   DMPlexGetStratumSize(dm, "vtk", 1, &numLabelCells);
268:   DMPlexGetStratumSize(dm, "vtk", 2, &numLabelVertices);
269:   hasLabel = numLabelCells > 0 || numLabelVertices > 0 ? PETSC_TRUE : PETSC_FALSE;
270:   for (p = pStart; p < pEnd; ++p) {
271:     /* Reject points not either cells or vertices */
272:     if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
273:     if (hasLabel) {
274:       PetscInt value;

276:       if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
277:           ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
278:         DMPlexGetLabelValue(dm, "vtk", p, &value);
279:         if (value != 1) continue;
280:       }
281:     }
282:     PetscSectionGetDof(section, p, &numDof);
283:     if (numDof) break;
284:   }
285:   MPI_Allreduce(&numDof, &maxDof, 1, MPIU_INT, MPI_MAX, comm);
286:   enforceDof = PetscMax(enforceDof, maxDof);
287:   VecGetArray(v, &array);
288:   if (!rank) {
289:     char formatString[8];

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

296:       /* Reject points not either cells or vertices */
297:       if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
298:       if (hasLabel) {
299:         PetscInt value;

301:         if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
302:             ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
303:           DMPlexGetLabelValue(dm, "vtk", p, &value);
304:           if (value != 1) continue;
305:         }
306:       }
307:       PetscSectionGetDof(section, p, &dof);
308:       PetscSectionGetOffset(section, p, &off);
309:       PetscSectionGetOffset(globalSection, p, &goff);
310:       if (dof && goff >= 0) {
311:         for (d = 0; d < dof; d++) {
312:           if (d > 0) {
313:             PetscFPrintf(comm, fp, " ");
314:           }
315:           PetscFPrintf(comm, fp, formatString, PetscRealPart(array[off+d])*scale);
316:         }
317:         for (d = dof; d < enforceDof; d++) {
318:           PetscFPrintf(comm, fp, " 0.0");
319:         }
320:         PetscFPrintf(comm, fp, "\n");
321:       }
322:     }
323:     for (proc = 1; proc < numProcs; ++proc) {
324:       PetscScalar *remoteValues;
325:       PetscInt    size = 0, d;
326:       MPI_Status  status;

328:       MPI_Recv(&size, 1, MPIU_INT, proc, tag, comm, &status);
329:       PetscMalloc(size * sizeof(PetscScalar), &remoteValues);
330:       MPI_Recv(remoteValues, size, mpiType, proc, tag, comm, &status);
331:       for (p = 0; p < size/maxDof; ++p) {
332:         for (d = 0; d < maxDof; ++d) {
333:           if (d > 0) {
334:             PetscFPrintf(comm, fp, " ");
335:           }
336:           PetscFPrintf(comm, fp, formatString, PetscRealPart(remoteValues[p*maxDof+d])*scale);
337:         }
338:         for (d = maxDof; d < enforceDof; ++d) {
339:           PetscFPrintf(comm, fp, " 0.0");
340:         }
341:         PetscFPrintf(comm, fp, "\n");
342:       }
343:       PetscFree(remoteValues);
344:     }
345:   } else {
346:     PetscScalar *localValues;
347:     PetscInt    size, k = 0;

349:     PetscSectionGetStorageSize(section, &size);
350:     PetscMalloc(size * sizeof(PetscScalar), &localValues);
351:     for (p = pStart; p < pEnd; ++p) {
352:       PetscInt dof, off, goff, d;

354:       /* Reject points not either cells or vertices */
355:       if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
356:       if (hasLabel) {
357:         PetscInt value;

359:         if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
360:             ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
361:           DMPlexGetLabelValue(dm, "vtk", p, &value);
362:           if (value != 1) continue;
363:         }
364:       }
365:       PetscSectionGetDof(section, p, &dof);
366:       PetscSectionGetOffset(section, p, &off);
367:       PetscSectionGetOffset(globalSection, p, &goff);
368:       if (goff >= 0) {
369:         for (d = 0; d < dof; ++d) {
370:           localValues[k++] = array[off+d];
371:         }
372:       }
373:     }
374:     MPI_Send(&k, 1, MPIU_INT, 0, tag, comm);
375:     MPI_Send(localValues, k, mpiType, 0, tag, comm);
376:     PetscFree(localValues);
377:   }
378:   VecRestoreArray(v, &array);
379:   return(0);
380: }

384: PetscErrorCode DMPlexVTKWriteField_ASCII(DM dm, PetscSection section, PetscSection globalSection, Vec field, const char name[], FILE *fp, PetscInt enforceDof, PetscInt precision, PetscReal scale)
385: {
386:   MPI_Comm       comm;
387:   PetscInt       numDof = 0, maxDof;
388:   PetscInt       pStart, pEnd, p;

392:   PetscObjectGetComm((PetscObject)dm,&comm);
393:   PetscSectionGetChart(section, &pStart, &pEnd);
394:   for (p = pStart; p < pEnd; ++p) {
395:     PetscSectionGetDof(section, p, &numDof);
396:     if (numDof) break;
397:   }
398:   numDof = PetscMax(numDof, enforceDof);
399:   MPI_Allreduce(&numDof, &maxDof, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
400:   if (!name) name = "Unknown";
401:   if (maxDof == 3) {
402:     PetscFPrintf(comm, fp, "VECTORS %s double\n", name);
403:   } else {
404:     PetscFPrintf(comm, fp, "SCALARS %s double %d\n", name, maxDof);
405:     PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n");
406:   }
407:   DMPlexVTKWriteSection_ASCII(dm, section, globalSection, field, fp, enforceDof, precision, scale);
408:   return(0);
409: }

413: static PetscErrorCode DMPlexVTKWriteAll_ASCII(DM dm, PetscViewer viewer)
414: {
415:   MPI_Comm                 comm;
416:   PetscViewer_VTK          *vtk = (PetscViewer_VTK*) viewer->data;
417:   FILE                     *fp;
418:   PetscViewerVTKObjectLink link;
419:   PetscSection             coordSection, globalCoordSection;
420:   PetscLayout              vLayout;
421:   Vec                      coordinates;
422:   PetscReal                lengthScale;
423:   PetscInt                 vMax, totVertices, totCells;
424:   PetscBool                hasPoint = PETSC_FALSE, hasCell = PETSC_FALSE;
425:   PetscErrorCode           ierr;

428:   PetscObjectGetComm((PetscObject)dm,&comm);
429:   PetscFOpen(comm, vtk->filename, "wb", &fp);
430:   PetscFPrintf(comm, fp, "# vtk DataFile Version 2.0\n");
431:   PetscFPrintf(comm, fp, "Simplicial Mesh Example\n");
432:   PetscFPrintf(comm, fp, "ASCII\n");
433:   PetscFPrintf(comm, fp, "DATASET UNSTRUCTURED_GRID\n");
434:   /* Vertices */
435:   DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
436:   DMPlexGetCoordinateSection(dm, &coordSection);
437:   PetscSectionCreateGlobalSection(coordSection, dm->sf, PETSC_FALSE, &globalCoordSection);
438:   DMGetCoordinatesLocal(dm, &coordinates);
439:   DMPlexGetHybridBounds(dm, NULL, NULL, NULL, &vMax);
440:   if (vMax >= 0) {
441:     PetscInt pStart, pEnd, p, localSize = 0;

443:     PetscSectionGetChart(globalCoordSection, &pStart, &pEnd);
444:     pEnd = PetscMin(pEnd, vMax);
445:     for (p = pStart; p < pEnd; ++p) {
446:       PetscInt dof;

448:       PetscSectionGetDof(globalCoordSection, p, &dof);
449:       if (dof > 0) ++localSize;
450:     }
451:     PetscLayoutCreate(PetscObjectComm((PetscObject)dm), &vLayout);
452:     PetscLayoutSetLocalSize(vLayout, localSize);
453:     PetscLayoutSetBlockSize(vLayout, 1);
454:     PetscLayoutSetUp(vLayout);
455:   } else {
456:     PetscSectionGetPointLayout(PetscObjectComm((PetscObject)dm), globalCoordSection, &vLayout);
457:   }
458:   PetscLayoutGetSize(vLayout, &totVertices);
459:   PetscFPrintf(comm, fp, "POINTS %d double\n", totVertices);
460:   DMPlexVTKWriteSection_ASCII(dm, coordSection, globalCoordSection, coordinates, fp, 3, PETSC_DETERMINE, lengthScale);
461:   /* Cells */
462:   DMPlexVTKWriteCells_ASCII(dm, fp, &totCells);
463:   /* Vertex fields */
464:   for (link = vtk->link; link; link = link->next) {
465:     if ((link->ft == PETSC_VTK_POINT_FIELD) || (link->ft == PETSC_VTK_POINT_VECTOR_FIELD)) hasPoint = PETSC_TRUE;
466:     if ((link->ft == PETSC_VTK_CELL_FIELD)  || (link->ft == PETSC_VTK_CELL_VECTOR_FIELD))  hasCell  = PETSC_TRUE;
467:   }
468:   if (hasPoint) {
469:     PetscFPrintf(comm, fp, "POINT_DATA %d\n", totVertices);
470:     for (link = vtk->link; link; link = link->next) {
471:       Vec          X = (Vec) link->vec;
472:       DM           dmX;
473:       PetscSection section, globalSection, newSection = NULL;
474:       const char   *name;
475:       PetscInt     enforceDof = PETSC_DETERMINE;

477:       if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
478:       if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) enforceDof = 3;
479:       PetscObjectGetName(link->vec, &name);
480:       VecGetDM(X, &dmX);
481:       if (dmX) {
482:         DMLabel  subpointMap, subpointMapX;
483:         PetscInt dim, dimX, pStart, pEnd, qStart, qEnd;

485:         DMGetDefaultSection(dmX, &section);
486:         /* Here is where we check whether dmX is a submesh of dm */
487:         DMPlexGetDimension(dm,  &dim);
488:         DMPlexGetDimension(dmX, &dimX);
489:         DMPlexGetChart(dm,  &pStart, &pEnd);
490:         DMPlexGetChart(dmX, &qStart, &qEnd);
491:         DMPlexGetSubpointMap(dm,  &subpointMap);
492:         DMPlexGetSubpointMap(dmX, &subpointMapX);
493:         if (((dim != dimX) || ((pEnd-pStart) < (qEnd-qStart))) && subpointMap && !subpointMapX) {
494:           const PetscInt *ind;
495:           IS              subpointIS;
496:           PetscInt        n, q;

498:           PetscPrintf(PETSC_COMM_SELF, "Making translation PetscSection\n");
499:           PetscSectionGetChart(section, &qStart, &qEnd);
500:           DMPlexCreateSubpointIS(dm, &subpointIS);
501:           ISGetLocalSize(subpointIS, &n);
502:           ISGetIndices(subpointIS, &ind);
503:           PetscSectionCreate(comm, &newSection);
504:           PetscSectionSetChart(newSection, pStart, pEnd);
505:           for (q = qStart; q < qEnd; ++q) {
506:             PetscInt dof, off, p;

508:             PetscSectionGetDof(section, q, &dof);
509:             if (dof) {
510:               PetscFindInt(q, n, ind, &p);
511:               if (p >= pStart) {
512:                 PetscSectionSetDof(newSection, p, dof);
513:                 PetscSectionGetOffset(section, q, &off);
514:                 PetscSectionSetOffset(newSection, p, off);
515:               }
516:             }
517:           }
518:           ISRestoreIndices(subpointIS, &ind);
519:           ISDestroy(&subpointIS);
520:           /* No need to setup section */
521:           section = newSection;
522:         }
523:       } else {
524:         PetscObjectQuery(link->vec, "section", (PetscObject*) &section);
525:         if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
526:       }
527:       if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
528:       PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, &globalSection);
529:       DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, name, fp, enforceDof, PETSC_DETERMINE, 1.0);
530:       PetscSectionDestroy(&globalSection);
531:       if (newSection) {PetscSectionDestroy(&newSection);}
532:     }
533:   }
534:   /* Cell Fields */
535:   if (hasCell) {
536:     PetscFPrintf(comm, fp, "CELL_DATA %d\n", totCells);
537:     for (link = vtk->link; link; link = link->next) {
538:       Vec          X = (Vec) link->vec;
539:       DM           dmX;
540:       PetscSection section, globalSection;
541:       const char   *name;
542:       PetscInt     enforceDof = PETSC_DETERMINE;

544:       if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
545:       if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) enforceDof = 3;
546:       PetscObjectGetName(link->vec, &name);
547:       VecGetDM(X, &dmX);
548:       if (dmX) {
549:         DMGetDefaultSection(dmX, &section);
550:       } else {
551:         PetscContainer c;

553:         PetscObjectQuery(link->vec, "section", (PetscObject*) &c);
554:         if (!c) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
555:         PetscContainerGetPointer(c, (void**) &section);
556:       }
557:       if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
558:       PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, &globalSection);
559:       DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, name, fp, enforceDof, PETSC_DETERMINE, 1.0);
560:       PetscSectionDestroy(&globalSection);
561:     }
562:   }
563:   /* Cleanup */
564:   PetscSectionDestroy(&globalCoordSection);
565:   PetscLayoutDestroy(&vLayout);
566:   PetscFClose(comm, fp);
567:   return(0);
568: }

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

575:   Collective

577:   Input Arguments:
578: + odm - The DMPlex specifying the mesh, passed as a PetscObject
579: - viewer - viewer of type VTK

581:   Level: developer

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

588: .seealso: PETSCVIEWERVTK
589: @*/
590: PetscErrorCode DMPlexVTKWriteAll(PetscObject odm, PetscViewer viewer)
591: {
592:   DM             dm = (DM) odm;
593:   PetscBool      isvtk;

599:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);
600:   if (!isvtk) SETERRQ1(PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_INCOMP, "Cannot use viewer type %s", ((PetscObject)viewer)->type_name);
601:   switch (viewer->format) {
602:   case PETSC_VIEWER_ASCII_VTK:
603:     DMPlexVTKWriteAll_ASCII(dm, viewer);
604:     break;
605:   case PETSC_VIEWER_VTK_VTU:
606:     DMPlexVTKWriteAll_VTU(dm, viewer);
607:     break;
608:   default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
609:   }
610:   return(0);
611: }