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:   PetscFunctionBegin;
      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 5:
     55:       *cellType = 14; /* VTK_PYRAMID */
     56:       break;
     57:     case 6:
     58:       *cellType = 13; /* VTK_WEDGE */
     59:       break;
     60:     case 8:
     61:       *cellType = 12; /* VTK_HEXAHEDRON */
     62:       break;
     63:     case 10:
     64:       *cellType = 24; /* VTK_QUADRATIC_TETRA */
     65:       break;
     66:     case 27:
     67:       *cellType = 29; /* VTK_QUADRATIC_HEXAHEDRON */
     68:       break;
     69:     default:
     70:       break;
     71:     }
     72:   }
     73:   PetscFunctionReturn(PETSC_SUCCESS);
     74: }

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

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

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

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

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

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

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

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

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

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

    210:       PetscCallMPI(MPIU_Recv(&numCells, 1, MPIU_INT, proc, tag, comm, &status));
    211:       PetscCallMPI(MPIU_Recv(corners, numCells, MPIU_INT, proc, tag, comm, &status));
    212:       for (c = 0; c < numCells; ++c) {
    213:         PetscCall(DMPlexVTKGetCellType_Internal(dm, dim, corners[c], &cellType));
    214:         PetscCall(PetscFPrintf(comm, fp, "%" PetscInt_FMT "\n", cellType));
    215:       }
    216:     }
    217:   } else {
    218:     PetscCallMPI(MPIU_Send(&numCells, 1, MPIU_INT, 0, tag, comm));
    219:     PetscCallMPI(MPIU_Send(corners, numCells, MPIU_INT, 0, tag, comm));
    220:   }
    221:   PetscCall(PetscFree(corners));
    222:   *totalCells = totCells;
    223:   PetscFunctionReturn(PETSC_SUCCESS);
    224: }

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

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

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

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

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

    274: static PetscErrorCode DMPlexVTKWriteSection_ASCII(DM dm, PetscSection section, PetscSection globalSection, Vec v, FILE *fp, PetscInt enforceDof, PetscInt precision, PetscReal scale, PetscInt imag)
    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, vStart, vEnd, pStart, pEnd, p;
    281:   PetscMPIInt        size, rank, proc, tag;
    282:   PetscBool          hasLabel;

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

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

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

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

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

    358:       PetscCallMPI(MPI_Recv(&size, 1, MPIU_INT, proc, tag, comm, &status));
    359:       PetscCall(PetscMalloc1(size, &remoteValues));
    360:       PetscCallMPI(MPIU_Recv(remoteValues, size, mpiType, proc, tag, comm, &status));
    361:       for (p = 0; p < size / maxDof; ++p) {
    362:         for (d = 0; d < maxDof; ++d) {
    363:           if (d > 0) PetscCall(PetscFPrintf(comm, fp, " "));
    364:           val  = remoteValues[p * maxDof + d];
    365:           dval = (PetscVTKReal)((imag ? PetscImaginaryPart(val) : PetscRealPart(val)) * scale);
    366:           PetscCall(PetscFPrintf(comm, fp, formatString, dval));
    367:         }
    368:         for (d = maxDof; d < enforceDof; ++d) PetscCall(PetscFPrintf(comm, fp, " 0.0"));
    369:         PetscCall(PetscFPrintf(comm, fp, "\n"));
    370:       }
    371:       PetscCall(PetscFree(remoteValues));
    372:     }
    373:   } else {
    374:     PetscScalar *localValues;
    375:     PetscInt     size, k = 0;

    377:     PetscCall(PetscSectionGetStorageSize(section, &size));
    378:     PetscCall(PetscMalloc1(size, &localValues));
    379:     for (p = pStart; p < pEnd; ++p) {
    380:       PetscInt dof, off, goff, d;

    382:       /* Reject points not either cells or vertices */
    383:       if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
    384:       if (hasLabel) {
    385:         PetscInt value;

    387:         if (((p >= cStart) && (p < cEnd) && numLabelCells) || ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
    388:           PetscCall(DMGetLabelValue(dm, "vtk", p, &value));
    389:           if (value != 1) continue;
    390:         }
    391:       }
    392:       PetscCall(PetscSectionGetDof(section, p, &dof));
    393:       PetscCall(PetscSectionGetOffset(section, p, &off));
    394:       PetscCall(PetscSectionGetOffset(globalSection, p, &goff));
    395:       if (goff >= 0) {
    396:         for (d = 0; d < dof; ++d) localValues[k++] = array[off + d];
    397:       }
    398:     }
    399:     PetscCallMPI(MPIU_Send(&k, 1, MPIU_INT, 0, tag, comm));
    400:     PetscCallMPI(MPIU_Send(localValues, k, mpiType, 0, tag, comm));
    401:     PetscCall(PetscFree(localValues));
    402:   }
    403:   PetscCall(VecRestoreArray(v, &array));
    404:   PetscFunctionReturn(PETSC_SUCCESS);
    405: }

    407: 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)
    408: {
    409:   MPI_Comm comm;
    410:   PetscInt numDof = 0, maxDof;
    411:   PetscInt pStart, pEnd, p;

    413:   PetscFunctionBegin;
    414:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
    415:   PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
    416:   for (p = pStart; p < pEnd; ++p) {
    417:     PetscCall(PetscSectionGetDof(section, p, &numDof));
    418:     if (numDof) break;
    419:   }
    420:   numDof = PetscMax(numDof, enforceDof);
    421:   PetscCallMPI(MPIU_Allreduce(&numDof, &maxDof, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm)));
    422:   if (!name) name = "Unknown";
    423:   if (maxDof == 3) {
    424:     if (nameComplex) {
    425:       PetscCall(PetscFPrintf(comm, fp, "VECTORS %s.%s double\n", name, imag ? "Im" : "Re"));
    426:     } else {
    427:       PetscCall(PetscFPrintf(comm, fp, "VECTORS %s double\n", name));
    428:     }
    429:   } else {
    430:     if (nameComplex) {
    431:       PetscCall(PetscFPrintf(comm, fp, "SCALARS %s.%s double %" PetscInt_FMT "\n", name, imag ? "Im" : "Re", maxDof));
    432:     } else {
    433:       PetscCall(PetscFPrintf(comm, fp, "SCALARS %s double %" PetscInt_FMT "\n", name, maxDof));
    434:     }
    435:     PetscCall(PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n"));
    436:   }
    437:   PetscCall(DMPlexVTKWriteSection_ASCII(dm, section, globalSection, field, fp, enforceDof, precision, scale, imag));
    438:   PetscFunctionReturn(PETSC_SUCCESS);
    439: }

    441: static PetscErrorCode DMPlexVTKWriteAll_ASCII(DM dm, PetscViewer viewer)
    442: {
    443:   MPI_Comm                 comm;
    444:   PetscViewer_VTK         *vtk = (PetscViewer_VTK *)viewer->data;
    445:   FILE                    *fp;
    446:   PetscViewerVTKObjectLink link;
    447:   PetscInt                 totVertices, totCells = 0, loops_per_scalar, l;
    448:   PetscBool                hasPoint = PETSC_FALSE, hasCell = PETSC_FALSE, writePartition = PETSC_FALSE, localized, writeComplex;
    449:   const char              *dmname;

    451:   PetscFunctionBegin;
    452: #if defined(PETSC_USE_COMPLEX)
    453:   loops_per_scalar = 2;
    454:   writeComplex     = PETSC_TRUE;
    455: #else
    456:   loops_per_scalar = 1;
    457:   writeComplex     = PETSC_FALSE;
    458: #endif
    459:   PetscCall(DMGetCoordinatesLocalized(dm, &localized));
    460:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
    461:   PetscCheck(!localized, comm, PETSC_ERR_SUP, "VTK output with localized coordinates not yet supported");
    462:   PetscCall(PetscFOpen(comm, vtk->filename, "wb", &fp));
    463:   PetscCall(PetscObjectGetName((PetscObject)dm, &dmname));
    464:   PetscCall(PetscFPrintf(comm, fp, "# vtk DataFile Version 2.0\n"));
    465:   PetscCall(PetscFPrintf(comm, fp, "%s\n", dmname));
    466:   PetscCall(PetscFPrintf(comm, fp, "ASCII\n"));
    467:   PetscCall(PetscFPrintf(comm, fp, "DATASET UNSTRUCTURED_GRID\n"));
    468:   /* Vertices */
    469:   {
    470:     PetscSection section, coordSection, globalCoordSection;
    471:     Vec          coordinates;
    472:     PetscReal    lengthScale;
    473:     DMLabel      label;
    474:     IS           vStratumIS;
    475:     PetscLayout  vLayout;

    477:     PetscCall(DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale));
    478:     PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
    479:     PetscCall(DMPlexGetDepthLabel(dm, &label));
    480:     PetscCall(DMLabelGetStratumIS(label, 0, &vStratumIS));
    481:     PetscCall(DMGetCoordinateSection(dm, &section));                                   /* This section includes all points */
    482:     PetscCall(PetscSectionCreateSubdomainSection(section, vStratumIS, &coordSection)); /* This one includes just vertices */
    483:     PetscCall(PetscSectionCreateGlobalSection(coordSection, dm->sf, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, &globalCoordSection));
    484:     PetscCall(PetscSectionGetPointLayout(comm, globalCoordSection, &vLayout));
    485:     PetscCall(PetscLayoutGetSize(vLayout, &totVertices));
    486:     PetscCall(PetscFPrintf(comm, fp, "POINTS %" PetscInt_FMT " double\n", totVertices));
    487:     PetscCall(DMPlexVTKWriteSection_ASCII(dm, coordSection, globalCoordSection, coordinates, fp, 3, PETSC_DETERMINE, lengthScale, 0));
    488:     PetscCall(ISDestroy(&vStratumIS));
    489:     PetscCall(PetscLayoutDestroy(&vLayout));
    490:     PetscCall(PetscSectionDestroy(&coordSection));
    491:     PetscCall(PetscSectionDestroy(&globalCoordSection));
    492:   }
    493:   /* Cells */
    494:   PetscCall(DMPlexVTKWriteCells_ASCII(dm, fp, &totCells));
    495:   /* Vertex fields */
    496:   for (link = vtk->link; link; link = link->next) {
    497:     if ((link->ft == PETSC_VTK_POINT_FIELD) || (link->ft == PETSC_VTK_POINT_VECTOR_FIELD)) hasPoint = PETSC_TRUE;
    498:     if ((link->ft == PETSC_VTK_CELL_FIELD) || (link->ft == PETSC_VTK_CELL_VECTOR_FIELD)) hasCell = PETSC_TRUE;
    499:   }
    500:   if (hasPoint) {
    501:     PetscCall(PetscFPrintf(comm, fp, "POINT_DATA %" PetscInt_FMT "\n", totVertices));
    502:     for (link = vtk->link; link; link = link->next) {
    503:       Vec          X       = (Vec)link->vec;
    504:       PetscSection section = NULL, globalSection, newSection = NULL;
    505:       char         namebuf[256];
    506:       const char  *name;
    507:       PetscInt     enforceDof = PETSC_DETERMINE;

    509:       if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
    510:       if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) enforceDof = 3;
    511:       PetscCall(PetscObjectGetName(link->vec, &name));
    512:       PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)&section));
    513:       if (!section) {
    514:         DM dmX;

    516:         PetscCall(VecGetDM(X, &dmX));
    517:         if (dmX) {
    518:           DMLabel  subpointMap, subpointMapX;
    519:           PetscInt dim, dimX, pStart, pEnd, qStart, qEnd;

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

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

    545:               PetscCall(PetscSectionGetDof(section, q, &dof));
    546:               if (dof) {
    547:                 PetscCall(PetscFindInt(q, n, ind, &p));
    548:                 if (p >= pStart) {
    549:                   PetscCall(PetscSectionSetDof(newSection, p, dof));
    550:                   PetscCall(PetscSectionGetOffset(section, q, &off));
    551:                   PetscCall(PetscSectionSetOffset(newSection, p, off));
    552:                 }
    553:               }
    554:             }
    555:             if (subpointIS) PetscCall(ISRestoreIndices(subpointIS, &ind));
    556:             /* No need to setup section */
    557:             section = newSection;
    558:           }
    559:         }
    560:       }
    561:       PetscCheck(section, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it and could not create one from VecGetDM()", name);
    562:       if (link->field >= 0) {
    563:         const char *fieldname;

    565:         PetscCall(PetscSectionGetFieldName(section, link->field, &fieldname));
    566:         PetscCall(PetscSectionGetField(section, link->field, &section));
    567:         if (fieldname) {
    568:           PetscCall(PetscSNPrintf(namebuf, sizeof(namebuf), "%s%s", name, fieldname));
    569:         } else {
    570:           PetscCall(PetscSNPrintf(namebuf, sizeof(namebuf), "%s%" PetscInt_FMT, name, link->field));
    571:         }
    572:       } else {
    573:         PetscCall(PetscSNPrintf(namebuf, sizeof(namebuf), "%s", name));
    574:       }
    575:       PetscCall(PetscViewerVTKSanitizeName_Internal(namebuf, sizeof(namebuf)));
    576:       PetscCall(PetscSectionCreateGlobalSection(section, dm->sf, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, &globalSection));
    577:       for (l = 0; l < loops_per_scalar; l++) PetscCall(DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, namebuf, fp, enforceDof, PETSC_DETERMINE, 1.0, writeComplex, l));
    578:       PetscCall(PetscSectionDestroy(&globalSection));
    579:       if (newSection) PetscCall(PetscSectionDestroy(&newSection));
    580:     }
    581:   }
    582:   /* Cell Fields */
    583:   PetscCall(PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dm_view_partition", &writePartition, NULL));
    584:   if (hasCell || writePartition) {
    585:     PetscCall(PetscFPrintf(comm, fp, "CELL_DATA %" PetscInt_FMT "\n", totCells));
    586:     for (link = vtk->link; link; link = link->next) {
    587:       Vec          X       = (Vec)link->vec;
    588:       PetscSection section = NULL, globalSection;
    589:       const char  *name    = "";
    590:       char         namebuf[256];
    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:       PetscCall(PetscObjectGetName(link->vec, &name));
    596:       PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)&section));
    597:       if (!section) {
    598:         DM dmX;

    600:         PetscCall(VecGetDM(X, &dmX));
    601:         if (dmX) PetscCall(DMGetLocalSection(dmX, &section));
    602:       }
    603:       PetscCheck(section, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it and could not create one from VecGetDM()", name);
    604:       if (link->field >= 0) {
    605:         const char *fieldname;

    607:         PetscCall(PetscSectionGetFieldName(section, link->field, &fieldname));
    608:         PetscCall(PetscSectionGetField(section, link->field, &section));
    609:         if (fieldname) {
    610:           PetscCall(PetscSNPrintf(namebuf, sizeof(namebuf), "%s%s", name, fieldname));
    611:         } else {
    612:           PetscCall(PetscSNPrintf(namebuf, sizeof(namebuf), "%s%" PetscInt_FMT, name, link->field));
    613:         }
    614:       } else {
    615:         PetscCall(PetscSNPrintf(namebuf, sizeof(namebuf), "%s", name));
    616:       }
    617:       PetscCall(PetscViewerVTKSanitizeName_Internal(namebuf, sizeof(namebuf)));
    618:       PetscCall(PetscSectionCreateGlobalSection(section, dm->sf, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, &globalSection));
    619:       for (l = 0; l < loops_per_scalar; l++) PetscCall(DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, namebuf, fp, enforceDof, PETSC_DETERMINE, 1.0, writeComplex, l));
    620:       PetscCall(PetscSectionDestroy(&globalSection));
    621:     }
    622:     if (writePartition) {
    623:       PetscCall(PetscFPrintf(comm, fp, "SCALARS partition int 1\n"));
    624:       PetscCall(PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n"));
    625:       PetscCall(DMPlexVTKWritePartition_ASCII(dm, fp));
    626:     }
    627:   }
    628:   /* Cleanup */
    629:   PetscCall(PetscFClose(comm, fp));
    630:   PetscFunctionReturn(PETSC_SUCCESS);
    631: }

    633: /*@
    634:   DMPlexVTKWriteAll - Write a file containing all the fields that have been provided to the viewer

    636:   Collective

    638:   Input Parameters:
    639: + odm    - The `DMPLEX` specifying the mesh, passed as a `PetscObject`
    640: - viewer - viewer of type `PETSCVIEWERVTK`

    642:   Level: developer

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

    649: .seealso: [](ch_unstructured), `DM`, `PETSCVIEWEREXODUSII`, `DMPLEX`, `PETSCVIEWERVTK`
    650: @*/
    651: PetscErrorCode DMPlexVTKWriteAll(PetscObject odm, PetscViewer viewer)
    652: {
    653:   DM dm = (DM)odm;

    655:   PetscFunctionBegin;
    658:   switch (viewer->format) {
    659:   case PETSC_VIEWER_ASCII_VTK_DEPRECATED:
    660:     PetscCall(DMPlexVTKWriteAll_ASCII(dm, viewer));
    661:     break;
    662:   case PETSC_VIEWER_VTK_VTU:
    663:     PetscCall(DMPlexVTKWriteAll_VTU(dm, viewer));
    664:     break;
    665:   default:
    666:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
    667:   }
    668:   PetscFunctionReturn(PETSC_SUCCESS);
    669: }