Actual source code: plexvtk.c

  1: #define PETSCDM_DLL
  2: #include <petsc/private/dmpleximpl.h>
  3: #include <../src/sys/classes/viewer/impls/vtk/vtkvimpl.h>

  5: PetscErrorCode DMPlexVTKGetCellType_Internal(DM dm, PetscInt dim, PetscInt corners, PetscInt *cellType)
  6: {
  7:   *cellType = -1;
  8:   switch (dim) {
  9:   case 0:
 10:     switch (corners) {
 11:     case 1:
 12:       *cellType = 1; /* VTK_VERTEX */
 13:       break;
 14:     default:
 15:       break;
 16:     }
 17:     break;
 18:   case 1:
 19:     switch (corners) {
 20:     case 2:
 21:       *cellType = 3; /* VTK_LINE */
 22:       break;
 23:     case 3:
 24:       *cellType = 21; /* VTK_QUADRATIC_EDGE */
 25:       break;
 26:     default:
 27:       break;
 28:     }
 29:     break;
 30:   case 2:
 31:     switch (corners) {
 32:     case 3:
 33:       *cellType = 5; /* VTK_TRIANGLE */
 34:       break;
 35:     case 4:
 36:       *cellType = 9; /* VTK_QUAD */
 37:       break;
 38:     case 6:
 39:       *cellType = 22; /* VTK_QUADRATIC_TRIANGLE */
 40:       break;
 41:     case 9:
 42:       *cellType = 23; /* VTK_QUADRATIC_QUAD */
 43:       break;
 44:     default:
 45:       break;
 46:     }
 47:     break;
 48:   case 3:
 49:     switch (corners) {
 50:     case 4:
 51:       *cellType = 10; /* VTK_TETRA */
 52:       break;
 53:     case 6:
 54:       *cellType = 13; /* VTK_WEDGE */
 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: }

 72: static PetscErrorCode DMPlexVTKWriteCells_ASCII(DM dm, FILE *fp, PetscInt *totalCells)
 73: {
 74:   MPI_Comm       comm;
 75:   DMLabel        label;
 76:   IS             globalVertexNumbers = NULL;
 77:   const PetscInt *gvertex;
 78:   PetscInt       dim;
 79:   PetscInt       numCorners = 0, totCorners = 0, maxCorners, *corners;
 80:   PetscInt       numCells   = 0, totCells   = 0, maxCells, cellHeight;
 81:   PetscInt       numLabelCells, maxLabelCells, cStart, cEnd, c, vStart, vEnd, v;
 82:   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:   DMGetLabel(dm, "vtk", &label);
 93:   DMGetStratumSize(dm, "vtk", 1, &numLabelCells);
 94:   MPIU_Allreduce(&numLabelCells, &maxLabelCells, 1, MPIU_INT, MPI_MAX, comm);
 95:   if (!maxLabelCells) label = NULL;
 96:   for (c = cStart; c < cEnd; ++c) {
 97:     PetscInt *closure = NULL;
 98:     PetscInt closureSize, value;

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

123:     PetscMalloc1(maxCorners, &vertices);
124:     for (c = cStart, numCells = 0; c < cEnd; ++c) {
125:       PetscInt *closure = NULL;
126:       PetscInt closureSize, value, nC = 0;

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

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

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

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

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

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

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

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

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

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

259:       MPI_Recv(&numCells, 1, MPIU_INT, proc, tag, comm, &status);
260:       for (c = 0; c < numCells; ++c) PetscFPrintf(comm, fp, "%d\n", proc);
261:     }
262:   } else {
263:     MPI_Send(&numCells, 1, MPIU_INT, 0, tag, comm);
264:   }
265:   return 0;
266: }

268: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
269: typedef double PetscVTKReal;
270: #elif defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL___FP16)
271: typedef float PetscVTKReal;
272: #else
273: typedef PetscReal PetscVTKReal;
274: #endif

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

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:   pStart   = PetscMax(PetscMin(cStart, vStart), pStart);
299:   pEnd     = PetscMin(PetscMax(cEnd,   vEnd),   pEnd);
300:   DMGetStratumSize(dm, "vtk", 1, &numLabelCells);
301:   DMGetStratumSize(dm, "vtk", 2, &numLabelVertices);
302:   hasLabel = numLabelCells > 0 || numLabelVertices > 0 ? PETSC_TRUE : PETSC_FALSE;
303:   for (p = pStart; p < pEnd; ++p) {
304:     /* Reject points not either cells or vertices */
305:     if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
306:     if (hasLabel) {
307:       PetscInt value;

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

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

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

336:         if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
337:             ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
338:           DMGetLabelValue(dm, "vtk", p, &value);
339:           if (value != 1) continue;
340:         }
341:       }
342:       PetscSectionGetDof(section, p, &dof);
343:       PetscSectionGetOffset(section, p, &off);
344:       PetscSectionGetOffset(globalSection, p, &goff);
345:       if (dof && goff >= 0) {
346:         for (d = 0; d < dof; d++) {
347:           if (d > 0) {
348:             PetscFPrintf(comm, fp, " ");
349:           }
350:           val = array[off+d];
351:           dval = (PetscVTKReal) ((imag ? PetscImaginaryPart(val) : PetscRealPart(val)) * scale);
352:           PetscFPrintf(comm, fp, formatString, dval);
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:           val = remoteValues[p*maxDof+d];
374:           dval = (PetscVTKReal) ((imag ? PetscImaginaryPart(val) : PetscRealPart(val)) * scale);
375:           PetscFPrintf(comm, fp, formatString, dval);
376:         }
377:         for (d = maxDof; d < enforceDof; ++d) {
378:           PetscFPrintf(comm, fp, " 0.0");
379:         }
380:         PetscFPrintf(comm, fp, "\n");
381:       }
382:       PetscFree(remoteValues);
383:     }
384:   } else {
385:     PetscScalar *localValues;
386:     PetscInt    size, k = 0;

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

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

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

421: static PetscErrorCode DMPlexVTKWriteField_ASCII(DM dm, PetscSection section, PetscSection globalSection, Vec field, const char name[], FILE *fp, PetscInt enforceDof, PetscInt precision, PetscReal scale, PetscBool nameComplex, PetscInt imag)
422: {
423:   MPI_Comm       comm;
424:   PetscInt       numDof = 0, maxDof;
425:   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:     if (nameComplex) {
438:       PetscFPrintf(comm, fp, "VECTORS %s.%s double\n", name, imag ? "Im" : "Re");
439:     } else {
440:       PetscFPrintf(comm, fp, "VECTORS %s double\n", name);
441:     }
442:   } else {
443:     if (nameComplex) {
444:       PetscFPrintf(comm, fp, "SCALARS %s.%s double %D\n", name, imag ? "Im" : "Re", maxDof);
445:     } else {
446:       PetscFPrintf(comm, fp, "SCALARS %s double %D\n", name, maxDof);
447:     }
448:     PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n");
449:   }
450:   DMPlexVTKWriteSection_ASCII(dm, section, globalSection, field, fp, enforceDof, precision, scale, imag);
451:   return 0;
452: }

454: static PetscErrorCode DMPlexVTKWriteAll_ASCII(DM dm, PetscViewer viewer)
455: {
456:   MPI_Comm                 comm;
457:   PetscViewer_VTK          *vtk = (PetscViewer_VTK*) viewer->data;
458:   FILE                     *fp;
459:   PetscViewerVTKObjectLink link;
460:   PetscInt                 totVertices, totCells = 0, loops_per_scalar, l;
461:   PetscBool                hasPoint = PETSC_FALSE, hasCell = PETSC_FALSE, writePartition = PETSC_FALSE, localized, writeComplex;
462:   const char               *dmname;

464: #if defined(PETSC_USE_COMPLEX)
465:   loops_per_scalar = 2;
466:   writeComplex = PETSC_TRUE;
467: #else
468:   loops_per_scalar = 1;
469:   writeComplex = PETSC_FALSE;
470: #endif
471:   DMGetCoordinatesLocalized(dm,&localized);
472:   PetscObjectGetComm((PetscObject)dm,&comm);
474:   PetscFOpen(comm, vtk->filename, "wb", &fp);
475:   PetscObjectGetName((PetscObject)dm, &dmname);
476:   PetscFPrintf(comm, fp, "# vtk DataFile Version 2.0\n");
477:   PetscFPrintf(comm, fp, "%s\n", dmname);
478:   PetscFPrintf(comm, fp, "ASCII\n");
479:   PetscFPrintf(comm, fp, "DATASET UNSTRUCTURED_GRID\n");
480:   /* Vertices */
481:   {
482:     PetscSection  section, coordSection, globalCoordSection;
483:     Vec           coordinates;
484:     PetscReal     lengthScale;
485:     DMLabel       label;
486:     IS            vStratumIS;
487:     PetscLayout   vLayout;

489:     DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
490:     DMGetCoordinatesLocal(dm, &coordinates);
491:     DMPlexGetDepthLabel(dm, &label);
492:     DMLabelGetStratumIS(label, 0, &vStratumIS);
493:     DMGetCoordinateSection(dm, &section);                                 /* This section includes all points */
494:     PetscSectionCreateSubmeshSection(section, vStratumIS, &coordSection); /* This one includes just vertices */
495:     PetscSectionCreateGlobalSection(coordSection, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalCoordSection);
496:     PetscSectionGetPointLayout(comm, globalCoordSection, &vLayout);
497:     PetscLayoutGetSize(vLayout, &totVertices);
498:     PetscFPrintf(comm, fp, "POINTS %D double\n", totVertices);
499:     DMPlexVTKWriteSection_ASCII(dm, coordSection, globalCoordSection, coordinates, fp, 3, PETSC_DETERMINE, lengthScale, 0);
500:     ISDestroy(&vStratumIS);
501:     PetscLayoutDestroy(&vLayout);
502:     PetscSectionDestroy(&coordSection);
503:     PetscSectionDestroy(&globalCoordSection);
504:   }
505:   /* Cells */
506:   DMPlexVTKWriteCells_ASCII(dm, fp, &totCells);
507:   /* Vertex fields */
508:   for (link = vtk->link; link; link = link->next) {
509:     if ((link->ft == PETSC_VTK_POINT_FIELD) || (link->ft == PETSC_VTK_POINT_VECTOR_FIELD)) hasPoint = PETSC_TRUE;
510:     if ((link->ft == PETSC_VTK_CELL_FIELD)  || (link->ft == PETSC_VTK_CELL_VECTOR_FIELD))  hasCell  = PETSC_TRUE;
511:   }
512:   if (hasPoint) {
513:     PetscFPrintf(comm, fp, "POINT_DATA %D\n", totVertices);
514:     for (link = vtk->link; link; link = link->next) {
515:       Vec          X = (Vec) link->vec;
516:       PetscSection section = NULL, globalSection, newSection = NULL;
517:       char         namebuf[256];
518:       const char   *name;
519:       PetscInt     enforceDof = PETSC_DETERMINE;

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

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

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

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

557:               PetscSectionGetDof(section, q, &dof);
558:               if (dof) {
559:                 PetscFindInt(q, n, ind, &p);
560:                 if (p >= pStart) {
561:                   PetscSectionSetDof(newSection, p, dof);
562:                   PetscSectionGetOffset(section, q, &off);
563:                   PetscSectionSetOffset(newSection, p, off);
564:                 }
565:               }
566:             }
567:             if (subpointIS) {
568:               ISRestoreIndices(subpointIS, &ind);
569:             }
570:             /* No need to setup section */
571:             section = newSection;
572:           }
573:         }
574:       }
576:       if (link->field >= 0) {
577:         const char *fieldname;

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

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

616:         VecGetDM(X, &dmX);
617:         if (dmX) {
618:           DMGetLocalSection(dmX, &section);
619:         }
620:       }
622:       if (link->field >= 0) {
623:         const char *fieldname;

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

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

656:   Collective

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

662:   Level: developer

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

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

678:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);
680:   switch (viewer->format) {
681:   case PETSC_VIEWER_ASCII_VTK_DEPRECATED:
682:     DMPlexVTKWriteAll_ASCII(dm, viewer);
683:     break;
684:   case PETSC_VIEWER_VTK_VTU:
685:     DMPlexVTKWriteAll_VTU(dm, viewer);
686:     break;
687:   default: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
688:   }
689:   return 0;
690: }