Actual source code: grvtk.c

      1: #include <petsc/private/dmdaimpl.h>
      2: /*
      3:    Note that the API for using PETSCVIEWERVTK is totally wrong since its use requires
      4:    including the private vtkvimpl.h file. The code should be refactored.
      5: */
      6: #include <../src/sys/classes/viewer/impls/vtk/vtkvimpl.h>

      8: /* Helper function which determines if any DMDA fields are named.  This is used
      9:    as a proxy for the user's intention to use DMDA fields as distinct
     10:    scalar-valued fields as opposed to a single vector-valued field */
     11: static PetscErrorCode DMDAGetFieldsNamed(DM da, PetscBool *fieldsnamed)
     12: {
     13:   PetscInt f, bs;

     15:   PetscFunctionBegin;
     16:   *fieldsnamed = PETSC_FALSE;
     17:   PetscCall(DMDAGetDof(da, &bs));
     18:   for (f = 0; f < bs; ++f) {
     19:     const char *fieldname;
     20:     PetscCall(DMDAGetFieldName(da, f, &fieldname));
     21:     if (fieldname) {
     22:       *fieldsnamed = PETSC_TRUE;
     23:       break;
     24:     }
     25:   }
     26:   PetscFunctionReturn(PETSC_SUCCESS);
     27: }

     29: static PetscErrorCode DMDAVTKWriteAll_VTS(DM da, PetscViewer viewer)
     30: {
     31:   const char *byte_order = PetscBinaryBigEndian() ? "BigEndian" : "LittleEndian";
     32: #if defined(PETSC_USE_REAL_SINGLE)
     33:   const char precision[] = "Float32";
     34: #elif defined(PETSC_USE_REAL_DOUBLE)
     35:   const char precision[] = "Float64";
     36: #else
     37:   const char precision[] = "UnknownPrecision";
     38: #endif
     39:   MPI_Comm                 comm;
     40:   Vec                      Coords;
     41:   PetscViewer_VTK         *vtk = (PetscViewer_VTK *)viewer->data;
     42:   PetscViewerVTKObjectLink link;
     43:   FILE                    *fp;
     44:   PetscMPIInt              rank, size, tag;
     45:   DMDALocalInfo            info;
     46:   PetscInt                 dim, mx, my, mz, cdim, bs, maxnnodes, maxbs, i, j, k;
     47:   PetscInt64               boffset;
     48:   PetscInt                 rloc[6], (*grloc)[6] = NULL;
     49:   PetscScalar             *array, *array2;

     51:   PetscFunctionBegin;
     52:   PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
     53:   PetscCheck(!PetscDefined(USE_COMPLEX), PETSC_COMM_SELF, PETSC_ERR_SUP, "Complex values not supported");
     54:   PetscCallMPI(MPI_Comm_size(comm, &size));
     55:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
     56:   PetscCall(DMDAGetInfo(da, &dim, &mx, &my, &mz, NULL, NULL, NULL, &bs, NULL, NULL, NULL, NULL, NULL));
     57:   PetscCall(DMDAGetLocalInfo(da, &info));
     58:   PetscCall(DMGetCoordinates(da, &Coords));
     59:   if (Coords) {
     60:     PetscInt csize;
     61:     PetscCall(VecGetSize(Coords, &csize));
     62:     PetscCheck(csize % (mx * my * mz) == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Coordinate vector size mismatch");
     63:     cdim = csize / (mx * my * mz);
     64:     PetscCheck(cdim >= dim && cdim <= 3, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Coordinate vector size mismatch");
     65:   } else {
     66:     cdim = dim;
     67:   }

     69:   PetscCall(PetscFOpen(comm, vtk->filename, "wb", &fp));
     70:   PetscCall(PetscFPrintf(comm, fp, "<?xml version=\"1.0\"?>\n"));
     71:   PetscCall(PetscFPrintf(comm, fp, "<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"%s\" header_type=\"UInt64\">\n", byte_order));
     72:   PetscCall(PetscFPrintf(comm, fp, "  <StructuredGrid WholeExtent=\"%d %" PetscInt_FMT " %d %" PetscInt_FMT " %d %" PetscInt_FMT "\">\n", 0, mx - 1, 0, my - 1, 0, mz - 1));

     74:   if (rank == 0) PetscCall(PetscMalloc1(size * 6, &grloc));
     75:   rloc[0] = info.xs;
     76:   rloc[1] = info.xm;
     77:   rloc[2] = info.ys;
     78:   rloc[3] = info.ym;
     79:   rloc[4] = info.zs;
     80:   rloc[5] = info.zm;
     81:   PetscCallMPI(MPI_Gather(rloc, 6, MPIU_INT, rank == 0 ? grloc[0] : NULL, 6, MPIU_INT, 0, comm));

     83:   /* Write XML header */
     84:   maxnnodes = 0; /* Used for the temporary array size on rank 0 */
     85:   maxbs     = 0; /* Used for the temporary array size on rank 0 */
     86:   boffset   = 0; /* Offset into binary file */
     87:   for (PetscMPIInt r = 0; r < size; r++) {
     88:     PetscInt xs = -1, xm = -1, ys = -1, ym = -1, zs = -1, zm = -1, nnodes = 0;

     90:     if (rank == 0) {
     91:       xs     = grloc[r][0];
     92:       xm     = grloc[r][1];
     93:       ys     = grloc[r][2];
     94:       ym     = grloc[r][3];
     95:       zs     = grloc[r][4];
     96:       zm     = grloc[r][5];
     97:       nnodes = xm * ym * zm;
     98:     }
     99:     maxnnodes = PetscMax(maxnnodes, nnodes);
    100:     PetscCall(PetscFPrintf(comm, fp, "    <Piece Extent=\"%" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\">\n", xs, xs + xm - 1, ys, ys + ym - 1, zs, zs + zm - 1));
    101:     PetscCall(PetscFPrintf(comm, fp, "      <Points>\n"));
    102:     PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"Position\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, boffset));
    103:     boffset += 3 * nnodes * sizeof(PetscScalar) + sizeof(PetscInt64);
    104:     PetscCall(PetscFPrintf(comm, fp, "      </Points>\n"));

    106:     PetscCall(PetscFPrintf(comm, fp, "      <PointData Scalars=\"ScalarPointData\">\n"));
    107:     for (link = vtk->link; link; link = link->next) {
    108:       Vec         X = (Vec)link->vec;
    109:       PetscInt    bs, f;
    110:       DM          daCurr;
    111:       PetscBool   fieldsnamed;
    112:       const char *vecname = "Unnamed";

    114:       PetscCall(VecGetDM(X, &daCurr));
    115:       PetscCall(DMDAGetInfo(daCurr, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &bs, NULL, NULL, NULL, NULL, NULL));
    116:       maxbs = PetscMax(maxbs, bs);

    118:       if (((PetscObject)X)->name || link != vtk->link) PetscCall(PetscObjectGetName((PetscObject)X, &vecname));

    120:       /* If any fields are named, add scalar fields. Otherwise, add a vector field */
    121:       PetscCall(DMDAGetFieldsNamed(daCurr, &fieldsnamed));
    122:       if (fieldsnamed) {
    123:         for (f = 0; f < bs; f++) {
    124:           char        buf[256];
    125:           const char *fieldname;
    126:           PetscCall(DMDAGetFieldName(daCurr, f, &fieldname));
    127:           if (!fieldname) {
    128:             PetscCall(PetscSNPrintf(buf, sizeof(buf), "%" PetscInt_FMT, f));
    129:             fieldname = buf;
    130:           }
    131:           PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s.%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset));
    132:           boffset += nnodes * sizeof(PetscScalar) + sizeof(PetscInt64);
    133:         }
    134:       } else {
    135:         PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s\" NumberOfComponents=\"%" PetscInt_FMT "\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, bs, boffset));
    136:         boffset += bs * nnodes * sizeof(PetscScalar) + sizeof(PetscInt64);
    137:       }
    138:     }
    139:     PetscCall(PetscFPrintf(comm, fp, "      </PointData>\n"));
    140:     PetscCall(PetscFPrintf(comm, fp, "    </Piece>\n"));
    141:   }
    142:   PetscCall(PetscFPrintf(comm, fp, "  </StructuredGrid>\n"));
    143:   PetscCall(PetscFPrintf(comm, fp, "  <AppendedData encoding=\"raw\">\n"));
    144:   PetscCall(PetscFPrintf(comm, fp, "_"));

    146:   /* Now write the arrays. */
    147:   tag = ((PetscObject)viewer)->tag;
    148:   PetscCall(PetscMalloc2(maxnnodes * PetscMax(3, maxbs), &array, maxnnodes * PetscMax(3, maxbs), &array2));
    149:   for (PetscMPIInt r = 0; r < size; r++) {
    150:     MPI_Status status;
    151:     PetscInt   xs = -1, xm = -1, ys = -1, ym = -1, zs = -1, zm = -1, nnodes = 0;

    153:     if (rank == 0) {
    154:       xs     = grloc[r][0];
    155:       xm     = grloc[r][1];
    156:       ys     = grloc[r][2];
    157:       ym     = grloc[r][3];
    158:       zs     = grloc[r][4];
    159:       zm     = grloc[r][5];
    160:       nnodes = xm * ym * zm;
    161:     } else if (r == rank) {
    162:       nnodes = info.xm * info.ym * info.zm;
    163:     }

    165:     /* Write the coordinates */
    166:     if (Coords) {
    167:       const PetscScalar *coords;

    169:       PetscCall(VecGetArrayRead(Coords, &coords));
    170:       if (rank == 0) {
    171:         if (r) {
    172:           PetscMPIInt nn;

    174:           PetscCallMPI(MPIU_Recv(array, nnodes * cdim, MPIU_SCALAR, r, tag, comm, &status));
    175:           PetscCallMPI(MPI_Get_count(&status, MPIU_SCALAR, &nn));
    176:           PetscCheck(nn == nnodes * cdim, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Array size mismatch");
    177:         } else PetscCall(PetscArraycpy(array, coords, nnodes * cdim));
    178:         /* Transpose coordinates to VTK (C-style) ordering */
    179:         for (k = 0; k < zm; k++) {
    180:           for (j = 0; j < ym; j++) {
    181:             for (i = 0; i < xm; i++) {
    182:               PetscInt Iloc        = i + xm * (j + ym * k);
    183:               array2[Iloc * 3 + 0] = array[Iloc * cdim + 0];
    184:               array2[Iloc * 3 + 1] = cdim > 1 ? array[Iloc * cdim + 1] : 0.0;
    185:               array2[Iloc * 3 + 2] = cdim > 2 ? array[Iloc * cdim + 2] : 0.0;
    186:             }
    187:           }
    188:         }
    189:       } else if (r == rank) {
    190:         PetscCallMPI(MPIU_Send((void *)coords, nnodes * cdim, MPIU_SCALAR, 0, tag, comm));
    191:       }
    192:       PetscCall(VecRestoreArrayRead(Coords, &coords));
    193:     } else { /* Fabricate some coordinates using grid index */
    194:       for (k = 0; k < zm; k++) {
    195:         for (j = 0; j < ym; j++) {
    196:           for (i = 0; i < xm; i++) {
    197:             PetscInt Iloc        = i + xm * (j + ym * k);
    198:             array2[Iloc * 3 + 0] = xs + i;
    199:             array2[Iloc * 3 + 1] = ys + j;
    200:             array2[Iloc * 3 + 2] = zs + k;
    201:           }
    202:         }
    203:       }
    204:     }
    205:     PetscCall(PetscViewerVTKFWrite(viewer, fp, array2, nnodes * 3, MPIU_SCALAR));

    207:     /* Write each of the objects queued up for this file */
    208:     for (link = vtk->link; link; link = link->next) {
    209:       Vec                X = (Vec)link->vec;
    210:       const PetscScalar *x;
    211:       PetscInt           bs, f;
    212:       DM                 daCurr;
    213:       PetscBool          fieldsnamed;
    214:       PetscCall(VecGetDM(X, &daCurr));
    215:       PetscCall(DMDAGetInfo(daCurr, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &bs, NULL, NULL, NULL, NULL, NULL));
    216:       PetscCall(VecGetArrayRead(X, &x));
    217:       if (rank == 0) {
    218:         if (r) {
    219:           PetscMPIInt nn;

    221:           PetscCallMPI(MPIU_Recv(array, nnodes * bs, MPIU_SCALAR, r, tag, comm, &status));
    222:           PetscCallMPI(MPI_Get_count(&status, MPIU_SCALAR, &nn));
    223:           PetscCheck(nn == nnodes * bs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Array size mismatch receiving from rank %d", r);
    224:         } else PetscCall(PetscArraycpy(array, x, nnodes * bs));

    226:         /* If any fields are named, add scalar fields. Otherwise, add a vector field */
    227:         PetscCall(DMDAGetFieldsNamed(daCurr, &fieldsnamed));
    228:         if (fieldsnamed) {
    229:           for (f = 0; f < bs; f++) {
    230:             /* Extract and transpose the f'th field */
    231:             for (k = 0; k < zm; k++) {
    232:               for (j = 0; j < ym; j++) {
    233:                 for (i = 0; i < xm; i++) {
    234:                   PetscInt Iloc = i + xm * (j + ym * k);
    235:                   array2[Iloc]  = array[Iloc * bs + f];
    236:                 }
    237:               }
    238:             }
    239:             PetscCall(PetscViewerVTKFWrite(viewer, fp, array2, nnodes, MPIU_SCALAR));
    240:           }
    241:         } else PetscCall(PetscViewerVTKFWrite(viewer, fp, array, bs * nnodes, MPIU_SCALAR));
    242:       } else if (r == rank) PetscCallMPI(MPIU_Send((void *)x, nnodes * bs, MPIU_SCALAR, 0, tag, comm));
    243:       PetscCall(VecRestoreArrayRead(X, &x));
    244:     }
    245:   }
    246:   PetscCall(PetscFree2(array, array2));
    247:   PetscCall(PetscFree(grloc));

    249:   PetscCall(PetscFPrintf(comm, fp, "\n </AppendedData>\n"));
    250:   PetscCall(PetscFPrintf(comm, fp, "</VTKFile>\n"));
    251:   PetscCall(PetscFClose(comm, fp));
    252:   PetscFunctionReturn(PETSC_SUCCESS);
    253: }

    255: static PetscErrorCode DMDAVTKWriteAll_VTR(DM da, PetscViewer viewer)
    256: {
    257:   const char *byte_order = PetscBinaryBigEndian() ? "BigEndian" : "LittleEndian";
    258: #if defined(PETSC_USE_REAL_SINGLE)
    259:   const char precision[] = "Float32";
    260: #elif defined(PETSC_USE_REAL_DOUBLE)
    261:   const char precision[] = "Float64";
    262: #else
    263:   const char precision[] = "UnknownPrecision";
    264: #endif
    265:   MPI_Comm                 comm;
    266:   PetscViewer_VTK         *vtk = (PetscViewer_VTK *)viewer->data;
    267:   PetscViewerVTKObjectLink link;
    268:   FILE                    *fp;
    269:   PetscMPIInt              rank, size, tag;
    270:   DMDALocalInfo            info;
    271:   PetscInt                 dim, mx, my, mz, maxnnodes, maxbs, i, j, k;
    272:   PetscInt64               boffset;
    273:   PetscInt                 rloc[6], (*grloc)[6] = NULL;
    274:   PetscScalar             *array, *array2;

    276:   PetscFunctionBegin;
    277:   PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
    278:   PetscCheck(!PetscDefined(USE_COMPLEX), PETSC_COMM_SELF, PETSC_ERR_SUP, "Complex values not supported");
    279:   PetscCallMPI(MPI_Comm_size(comm, &size));
    280:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    281:   PetscCall(DMDAGetInfo(da, &dim, &mx, &my, &mz, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
    282:   PetscCall(DMDAGetLocalInfo(da, &info));
    283:   PetscCall(PetscFOpen(comm, vtk->filename, "wb", &fp));
    284:   PetscCall(PetscFPrintf(comm, fp, "<?xml version=\"1.0\"?>\n"));
    285:   PetscCall(PetscFPrintf(comm, fp, "<VTKFile type=\"RectilinearGrid\" version=\"0.1\" byte_order=\"%s\" header_type=\"UInt64\">\n", byte_order));
    286:   PetscCall(PetscFPrintf(comm, fp, "  <RectilinearGrid WholeExtent=\"%d %" PetscInt_FMT " %d %" PetscInt_FMT " %d %" PetscInt_FMT "\">\n", 0, mx - 1, 0, my - 1, 0, mz - 1));

    288:   if (rank == 0) PetscCall(PetscMalloc1(size * 6, &grloc));
    289:   rloc[0] = info.xs;
    290:   rloc[1] = info.xm;
    291:   rloc[2] = info.ys;
    292:   rloc[3] = info.ym;
    293:   rloc[4] = info.zs;
    294:   rloc[5] = info.zm;
    295:   PetscCallMPI(MPI_Gather(rloc, 6, MPIU_INT, rank == 0 ? grloc[0] : NULL, 6, MPIU_INT, 0, comm));

    297:   /* Write XML header */
    298:   maxnnodes = 0; /* Used for the temporary array size on rank 0 */
    299:   maxbs     = 0; /* Used for the temporary array size on rank 0 */
    300:   boffset   = 0; /* Offset into binary file */
    301:   for (PetscMPIInt r = 0; r < size; r++) {
    302:     PetscInt xs = -1, xm = -1, ys = -1, ym = -1, zs = -1, zm = -1, nnodes = 0;

    304:     if (rank == 0) {
    305:       xs     = grloc[r][0];
    306:       xm     = grloc[r][1];
    307:       ys     = grloc[r][2];
    308:       ym     = grloc[r][3];
    309:       zs     = grloc[r][4];
    310:       zm     = grloc[r][5];
    311:       nnodes = xm * ym * zm;
    312:     }
    313:     maxnnodes = PetscMax(maxnnodes, nnodes);
    314:     PetscCall(PetscFPrintf(comm, fp, "    <Piece Extent=\"%" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\">\n", xs, xs + xm - 1, ys, ys + ym - 1, zs, zs + zm - 1));
    315:     PetscCall(PetscFPrintf(comm, fp, "      <Coordinates>\n"));
    316:     PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"Xcoord\"  format=\"appended\"  offset=\"%" PetscInt64_FMT "\" />\n", precision, boffset));
    317:     boffset += xm * sizeof(PetscScalar) + sizeof(PetscInt64);
    318:     PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"Ycoord\"  format=\"appended\"  offset=\"%" PetscInt64_FMT "\" />\n", precision, boffset));
    319:     boffset += ym * sizeof(PetscScalar) + sizeof(PetscInt64);
    320:     PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"Zcoord\"  format=\"appended\"  offset=\"%" PetscInt64_FMT "\" />\n", precision, boffset));
    321:     boffset += zm * sizeof(PetscScalar) + sizeof(PetscInt64);
    322:     PetscCall(PetscFPrintf(comm, fp, "      </Coordinates>\n"));
    323:     PetscCall(PetscFPrintf(comm, fp, "      <PointData Scalars=\"ScalarPointData\">\n"));
    324:     for (link = vtk->link; link; link = link->next) {
    325:       Vec         X = (Vec)link->vec;
    326:       PetscInt    bs, f;
    327:       DM          daCurr;
    328:       PetscBool   fieldsnamed;
    329:       const char *vecname = "Unnamed";

    331:       PetscCall(VecGetDM(X, &daCurr));
    332:       PetscCall(DMDAGetInfo(daCurr, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &bs, NULL, NULL, NULL, NULL, NULL));
    333:       maxbs = PetscMax(maxbs, bs);
    334:       if (((PetscObject)X)->name || link != vtk->link) PetscCall(PetscObjectGetName((PetscObject)X, &vecname));

    336:       /* If any fields are named, add scalar fields. Otherwise, add a vector field */
    337:       PetscCall(DMDAGetFieldsNamed(daCurr, &fieldsnamed));
    338:       if (fieldsnamed) {
    339:         for (f = 0; f < bs; f++) {
    340:           char        buf[256];
    341:           const char *fieldname;
    342:           PetscCall(DMDAGetFieldName(daCurr, f, &fieldname));
    343:           if (!fieldname) {
    344:             PetscCall(PetscSNPrintf(buf, sizeof(buf), "%" PetscInt_FMT, f));
    345:             fieldname = buf;
    346:           }
    347:           PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s.%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset));
    348:           boffset += nnodes * sizeof(PetscScalar) + sizeof(PetscInt64);
    349:         }
    350:       } else {
    351:         PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s\" NumberOfComponents=\"%" PetscInt_FMT "\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, bs, boffset));
    352:         boffset += bs * nnodes * sizeof(PetscScalar) + sizeof(PetscInt64);
    353:       }
    354:     }
    355:     PetscCall(PetscFPrintf(comm, fp, "      </PointData>\n"));
    356:     PetscCall(PetscFPrintf(comm, fp, "    </Piece>\n"));
    357:   }
    358:   PetscCall(PetscFPrintf(comm, fp, "  </RectilinearGrid>\n"));
    359:   PetscCall(PetscFPrintf(comm, fp, "  <AppendedData encoding=\"raw\">\n"));
    360:   PetscCall(PetscFPrintf(comm, fp, "_"));

    362:   /* Now write the arrays. */
    363:   tag = ((PetscObject)viewer)->tag;
    364:   PetscCall(PetscMalloc2(PetscMax(maxnnodes * maxbs, info.xm + info.ym + info.zm), &array, PetscMax(maxnnodes * maxbs, info.xm + info.ym + info.zm), &array2));
    365:   for (PetscMPIInt r = 0; r < size; r++) {
    366:     MPI_Status status;
    367:     PetscInt   xs = -1, xm = -1, ys = -1, ym = -1, zs = -1, zm = -1, nnodes = 0;

    369:     if (rank == 0) {
    370:       xs     = grloc[r][0];
    371:       xm     = grloc[r][1];
    372:       ys     = grloc[r][2];
    373:       ym     = grloc[r][3];
    374:       zs     = grloc[r][4];
    375:       zm     = grloc[r][5];
    376:       nnodes = xm * ym * zm;
    377:     } else if (r == rank) {
    378:       nnodes = info.xm * info.ym * info.zm;
    379:     }
    380:     { /* Write the coordinates */
    381:       Vec Coords;

    383:       PetscCall(DMGetCoordinates(da, &Coords));
    384:       if (Coords) {
    385:         const PetscScalar *coords;
    386:         PetscCall(VecGetArrayRead(Coords, &coords));
    387:         if (rank == 0) {
    388:           if (r) {
    389:             PetscMPIInt nn;

    391:             PetscCallMPI(MPIU_Recv(array, xm + ym + zm, MPIU_SCALAR, r, tag, comm, &status));
    392:             PetscCallMPI(MPI_Get_count(&status, MPIU_SCALAR, &nn));
    393:             PetscCheck(nn == xm + ym + zm, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Array size mismatch");
    394:           } else {
    395:             /* x coordinates */
    396:             for (j = 0, k = 0, i = 0; i < xm; i++) {
    397:               PetscInt Iloc = i + xm * (j + ym * k);

    399:               array[i] = coords[Iloc * dim + 0];
    400:             }
    401:             /* y coordinates */
    402:             for (i = 0, k = 0, j = 0; j < ym; j++) {
    403:               PetscInt Iloc = i + xm * (j + ym * k);

    405:               array[j + xm] = dim > 1 ? coords[Iloc * dim + 1] : 0;
    406:             }
    407:             /* z coordinates */
    408:             for (i = 0, j = 0, k = 0; k < zm; k++) {
    409:               PetscInt Iloc = i + xm * (j + ym * k);

    411:               array[k + xm + ym] = dim > 2 ? coords[Iloc * dim + 2] : 0;
    412:             }
    413:           }
    414:         } else if (r == rank) {
    415:           xm = info.xm;
    416:           ym = info.ym;
    417:           zm = info.zm;
    418:           for (j = 0, k = 0, i = 0; i < xm; i++) {
    419:             PetscInt Iloc = i + xm * (j + ym * k);

    421:             array2[i] = coords[Iloc * dim + 0];
    422:           }
    423:           for (i = 0, k = 0, j = 0; j < ym; j++) {
    424:             PetscInt Iloc = i + xm * (j + ym * k);

    426:             array2[j + xm] = dim > 1 ? coords[Iloc * dim + 1] : 0;
    427:           }
    428:           for (i = 0, j = 0, k = 0; k < zm; k++) {
    429:             PetscInt Iloc = i + xm * (j + ym * k);

    431:             array2[k + xm + ym] = dim > 2 ? coords[Iloc * dim + 2] : 0;
    432:           }
    433:           PetscCallMPI(MPIU_Send((void *)array2, xm + ym + zm, MPIU_SCALAR, 0, tag, comm));
    434:         }
    435:         PetscCall(VecRestoreArrayRead(Coords, &coords));
    436:       } else { /* Fabricate some coordinates using grid index */
    437:         for (i = 0; i < xm; i++) array[i] = xs + i;
    438:         for (j = 0; j < ym; j++) array[j + xm] = ys + j;
    439:         for (k = 0; k < zm; k++) array[k + xm + ym] = zs + k;
    440:       }
    441:       if (rank == 0) {
    442:         PetscCall(PetscViewerVTKFWrite(viewer, fp, &array[0], xm, MPIU_SCALAR));
    443:         PetscCall(PetscViewerVTKFWrite(viewer, fp, &array[xm], ym, MPIU_SCALAR));
    444:         PetscCall(PetscViewerVTKFWrite(viewer, fp, &array[xm + ym], zm, MPIU_SCALAR));
    445:       }
    446:     }

    448:     /* Write each of the objects queued up for this file */
    449:     for (link = vtk->link; link; link = link->next) {
    450:       Vec                X = (Vec)link->vec;
    451:       const PetscScalar *x;
    452:       PetscInt           bs, f;
    453:       DM                 daCurr;
    454:       PetscBool          fieldsnamed;

    456:       PetscCall(VecGetDM(X, &daCurr));
    457:       PetscCall(DMDAGetInfo(daCurr, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &bs, NULL, NULL, NULL, NULL, NULL));

    459:       PetscCall(VecGetArrayRead(X, &x));
    460:       if (rank == 0) {
    461:         if (r) {
    462:           PetscMPIInt nn;

    464:           PetscCallMPI(MPIU_Recv(array, nnodes * bs, MPIU_SCALAR, r, tag, comm, &status));
    465:           PetscCallMPI(MPI_Get_count(&status, MPIU_SCALAR, &nn));
    466:           PetscCheck(nn == nnodes * bs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Array size mismatch receiving from rank %d", r);
    467:         } else PetscCall(PetscArraycpy(array, x, nnodes * bs));
    468:         /* If any fields are named, add scalar fields. Otherwise, add a vector field */
    469:         PetscCall(DMDAGetFieldsNamed(daCurr, &fieldsnamed));
    470:         if (fieldsnamed) {
    471:           for (f = 0; f < bs; f++) {
    472:             /* Extract and transpose the f'th field */
    473:             for (k = 0; k < zm; k++) {
    474:               for (j = 0; j < ym; j++) {
    475:                 for (i = 0; i < xm; i++) {
    476:                   PetscInt Iloc = i + xm * (j + ym * k);
    477:                   array2[Iloc]  = array[Iloc * bs + f];
    478:                 }
    479:               }
    480:             }
    481:             PetscCall(PetscViewerVTKFWrite(viewer, fp, array2, nnodes, MPIU_SCALAR));
    482:           }
    483:         }
    484:         PetscCall(PetscViewerVTKFWrite(viewer, fp, array, nnodes * bs, MPIU_SCALAR));

    486:       } else if (r == rank) {
    487:         PetscCallMPI(MPIU_Send((void *)x, nnodes * bs, MPIU_SCALAR, 0, tag, comm));
    488:       }
    489:       PetscCall(VecRestoreArrayRead(X, &x));
    490:     }
    491:   }
    492:   PetscCall(PetscFree2(array, array2));
    493:   PetscCall(PetscFree(grloc));

    495:   PetscCall(PetscFPrintf(comm, fp, "\n </AppendedData>\n"));
    496:   PetscCall(PetscFPrintf(comm, fp, "</VTKFile>\n"));
    497:   PetscCall(PetscFClose(comm, fp));
    498:   PetscFunctionReturn(PETSC_SUCCESS);
    499: }

    501: /*@
    502:   DMDAVTKWriteAll - Write a file containing all the fields that have been provided to the viewer

    504:   Collective

    506:   Input Parameters:
    507: + odm    - `DMDA` specifying the grid layout, passed as a `PetscObject`
    508: - viewer - viewer of type `PETSCVIEWERVTK`

    510:   Level: developer

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

    517:   If any fields have been named (see e.g. `DMDASetFieldName()`), then individual scalar
    518:   fields are written. Otherwise, a single multi-dof (vector) field is written.

    520: .seealso: [](sec_struct), `DMDA`, `DM`, `PETSCVIEWERVTK`, `DMDASetFieldName()`
    521: @*/
    522: PetscErrorCode DMDAVTKWriteAll(PetscObject odm, PetscViewer viewer)
    523: {
    524:   DM        dm = (DM)odm;
    525:   PetscBool isvtk;

    527:   PetscFunctionBegin;
    530:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk));
    531:   PetscCheck(isvtk, PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_INCOMP, "Cannot use viewer type %s", ((PetscObject)viewer)->type_name);
    532:   switch (viewer->format) {
    533:   case PETSC_VIEWER_VTK_VTS:
    534:     PetscCall(DMDAVTKWriteAll_VTS(dm, viewer));
    535:     break;
    536:   case PETSC_VIEWER_VTK_VTR:
    537:     PetscCall(DMDAVTKWriteAll_VTR(dm, viewer));
    538:     break;
    539:   default:
    540:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
    541:   }
    542:   PetscFunctionReturn(PETSC_SUCCESS);
    543: }