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

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

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

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

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

145:   /* Now write the arrays. */
146:   tag = ((PetscObject)viewer)->tag;
147:   PetscCall(PetscMalloc2(maxnnodes * PetscMax(3, maxbs), &array, maxnnodes * PetscMax(3, maxbs), &array2));
148:   for (r = 0; r < size; r++) {
149:     MPI_Status status;
150:     PetscInt   xs = -1, xm = -1, ys = -1, ym = -1, zs = -1, zm = -1, nnodes = 0;
151:     if (rank == 0) {
152:       xs     = grloc[r][0];
153:       xm     = grloc[r][1];
154:       ys     = grloc[r][2];
155:       ym     = grloc[r][3];
156:       zs     = grloc[r][4];
157:       zm     = grloc[r][5];
158:       nnodes = xm * ym * zm;
159:     } else if (r == rank) {
160:       nnodes = info.xm * info.ym * info.zm;
161:     }

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

203:     /* Write each of the objects queued up for this file */
204:     for (link = vtk->link; link; link = link->next) {
205:       Vec                X = (Vec)link->vec;
206:       const PetscScalar *x;
207:       PetscInt           bs, f;
208:       DM                 daCurr;
209:       PetscBool          fieldsnamed;
210:       PetscCall(VecGetDM(X, &daCurr));
211:       PetscCall(DMDAGetInfo(daCurr, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &bs, NULL, NULL, NULL, NULL, NULL));
212:       PetscCall(VecGetArrayRead(X, &x));
213:       if (rank == 0) {
214:         if (r) {
215:           PetscMPIInt nn;
216:           PetscCallMPI(MPI_Recv(array, nnodes * bs, MPIU_SCALAR, r, tag, comm, &status));
217:           PetscCallMPI(MPI_Get_count(&status, MPIU_SCALAR, &nn));
218:           PetscCheck(nn == nnodes * bs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Array size mismatch receiving from rank %" PetscInt_FMT, r);
219:         } else PetscCall(PetscArraycpy(array, x, nnodes * bs));

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

244:   PetscCall(PetscFPrintf(comm, fp, "\n </AppendedData>\n"));
245:   PetscCall(PetscFPrintf(comm, fp, "</VTKFile>\n"));
246:   PetscCall(PetscFClose(comm, fp));
247:   PetscFunctionReturn(PETSC_SUCCESS);
248: }

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

271:   PetscFunctionBegin;
272:   PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
273:   PetscCheck(!PetscDefined(USE_COMPLEX), PETSC_COMM_SELF, PETSC_ERR_SUP, "Complex values not supported");
274:   PetscCallMPI(MPI_Comm_size(comm, &size));
275:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
276:   PetscCall(DMDAGetInfo(da, &dim, &mx, &my, &mz, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
277:   PetscCall(DMDAGetLocalInfo(da, &info));
278:   PetscCall(PetscFOpen(comm, vtk->filename, "wb", &fp));
279:   PetscCall(PetscFPrintf(comm, fp, "<?xml version=\"1.0\"?>\n"));
280:   PetscCall(PetscFPrintf(comm, fp, "<VTKFile type=\"RectilinearGrid\" version=\"0.1\" byte_order=\"%s\" header_type=\"UInt64\">\n", byte_order));
281:   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));

283:   if (rank == 0) PetscCall(PetscMalloc1(size * 6, &grloc));
284:   rloc[0] = info.xs;
285:   rloc[1] = info.xm;
286:   rloc[2] = info.ys;
287:   rloc[3] = info.ym;
288:   rloc[4] = info.zs;
289:   rloc[5] = info.zm;
290:   PetscCallMPI(MPI_Gather(rloc, 6, MPIU_INT, &grloc[0][0], 6, MPIU_INT, 0, comm));

292:   /* Write XML header */
293:   maxnnodes = 0; /* Used for the temporary array size on rank 0 */
294:   maxbs     = 0; /* Used for the temporary array size on rank 0 */
295:   boffset   = 0; /* Offset into binary file */
296:   for (r = 0; r < size; r++) {
297:     PetscInt xs = -1, xm = -1, ys = -1, ym = -1, zs = -1, zm = -1, nnodes = 0;
298:     if (rank == 0) {
299:       xs     = grloc[r][0];
300:       xm     = grloc[r][1];
301:       ys     = grloc[r][2];
302:       ym     = grloc[r][3];
303:       zs     = grloc[r][4];
304:       zm     = grloc[r][5];
305:       nnodes = xm * ym * zm;
306:     }
307:     maxnnodes = PetscMax(maxnnodes, nnodes);
308:     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));
309:     PetscCall(PetscFPrintf(comm, fp, "      <Coordinates>\n"));
310:     PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"Xcoord\"  format=\"appended\"  offset=\"%" PetscInt64_FMT "\" />\n", precision, boffset));
311:     boffset += xm * sizeof(PetscScalar) + sizeof(PetscInt64);
312:     PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"Ycoord\"  format=\"appended\"  offset=\"%" PetscInt64_FMT "\" />\n", precision, boffset));
313:     boffset += ym * sizeof(PetscScalar) + sizeof(PetscInt64);
314:     PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"Zcoord\"  format=\"appended\"  offset=\"%" PetscInt64_FMT "\" />\n", precision, boffset));
315:     boffset += zm * sizeof(PetscScalar) + sizeof(PetscInt64);
316:     PetscCall(PetscFPrintf(comm, fp, "      </Coordinates>\n"));
317:     PetscCall(PetscFPrintf(comm, fp, "      <PointData Scalars=\"ScalarPointData\">\n"));
318:     for (link = vtk->link; link; link = link->next) {
319:       Vec         X = (Vec)link->vec;
320:       PetscInt    bs, f;
321:       DM          daCurr;
322:       PetscBool   fieldsnamed;
323:       const char *vecname = "Unnamed";

325:       PetscCall(VecGetDM(X, &daCurr));
326:       PetscCall(DMDAGetInfo(daCurr, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &bs, NULL, NULL, NULL, NULL, NULL));
327:       maxbs = PetscMax(maxbs, bs);
328:       if (((PetscObject)X)->name || link != vtk->link) PetscCall(PetscObjectGetName((PetscObject)X, &vecname));

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

356:   /* Now write the arrays. */
357:   tag = ((PetscObject)viewer)->tag;
358:   PetscCall(PetscMalloc2(PetscMax(maxnnodes * maxbs, info.xm + info.ym + info.zm), &array, PetscMax(maxnnodes * maxbs, info.xm + info.ym + info.zm), &array2));
359:   for (r = 0; r < size; r++) {
360:     MPI_Status status;
361:     PetscInt   xs = -1, xm = -1, ys = -1, ym = -1, zs = -1, zm = -1, nnodes = 0;
362:     if (rank == 0) {
363:       xs     = grloc[r][0];
364:       xm     = grloc[r][1];
365:       ys     = grloc[r][2];
366:       ym     = grloc[r][3];
367:       zs     = grloc[r][4];
368:       zm     = grloc[r][5];
369:       nnodes = xm * ym * zm;
370:     } else if (r == rank) {
371:       nnodes = info.xm * info.ym * info.zm;
372:     }
373:     { /* Write the coordinates */
374:       Vec Coords;

376:       PetscCall(DMGetCoordinates(da, &Coords));
377:       if (Coords) {
378:         const PetscScalar *coords;
379:         PetscCall(VecGetArrayRead(Coords, &coords));
380:         if (rank == 0) {
381:           if (r) {
382:             PetscMPIInt nn;
383:             PetscCallMPI(MPI_Recv(array, xm + ym + zm, MPIU_SCALAR, r, tag, comm, &status));
384:             PetscCallMPI(MPI_Get_count(&status, MPIU_SCALAR, &nn));
385:             PetscCheck(nn == xm + ym + zm, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Array size mismatch");
386:           } else {
387:             /* x coordinates */
388:             for (j = 0, k = 0, i = 0; i < xm; i++) {
389:               PetscInt Iloc = i + xm * (j + ym * k);
390:               array[i]      = coords[Iloc * dim + 0];
391:             }
392:             /* y coordinates */
393:             for (i = 0, k = 0, j = 0; j < ym; j++) {
394:               PetscInt Iloc = i + xm * (j + ym * k);
395:               array[j + xm] = dim > 1 ? coords[Iloc * dim + 1] : 0;
396:             }
397:             /* z coordinates */
398:             for (i = 0, j = 0, k = 0; k < zm; k++) {
399:               PetscInt Iloc      = i + xm * (j + ym * k);
400:               array[k + xm + ym] = dim > 2 ? coords[Iloc * dim + 2] : 0;
401:             }
402:           }
403:         } else if (r == rank) {
404:           xm = info.xm;
405:           ym = info.ym;
406:           zm = info.zm;
407:           for (j = 0, k = 0, i = 0; i < xm; i++) {
408:             PetscInt Iloc = i + xm * (j + ym * k);
409:             array2[i]     = coords[Iloc * dim + 0];
410:           }
411:           for (i = 0, k = 0, j = 0; j < ym; j++) {
412:             PetscInt Iloc  = i + xm * (j + ym * k);
413:             array2[j + xm] = dim > 1 ? coords[Iloc * dim + 1] : 0;
414:           }
415:           for (i = 0, j = 0, k = 0; k < zm; k++) {
416:             PetscInt Iloc       = i + xm * (j + ym * k);
417:             array2[k + xm + ym] = dim > 2 ? coords[Iloc * dim + 2] : 0;
418:           }
419:           PetscCallMPI(MPI_Send((void *)array2, xm + ym + zm, MPIU_SCALAR, 0, tag, comm));
420:         }
421:         PetscCall(VecRestoreArrayRead(Coords, &coords));
422:       } else { /* Fabricate some coordinates using grid index */
423:         for (i = 0; i < xm; i++) array[i] = xs + i;
424:         for (j = 0; j < ym; j++) array[j + xm] = ys + j;
425:         for (k = 0; k < zm; k++) array[k + xm + ym] = zs + k;
426:       }
427:       if (rank == 0) {
428:         PetscCall(PetscViewerVTKFWrite(viewer, fp, &(array[0]), xm, MPIU_SCALAR));
429:         PetscCall(PetscViewerVTKFWrite(viewer, fp, &(array[xm]), ym, MPIU_SCALAR));
430:         PetscCall(PetscViewerVTKFWrite(viewer, fp, &(array[xm + ym]), zm, MPIU_SCALAR));
431:       }
432:     }

434:     /* Write each of the objects queued up for this file */
435:     for (link = vtk->link; link; link = link->next) {
436:       Vec                X = (Vec)link->vec;
437:       const PetscScalar *x;
438:       PetscInt           bs, f;
439:       DM                 daCurr;
440:       PetscBool          fieldsnamed;
441:       PetscCall(VecGetDM(X, &daCurr));
442:       PetscCall(DMDAGetInfo(daCurr, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &bs, NULL, NULL, NULL, NULL, NULL));

444:       PetscCall(VecGetArrayRead(X, &x));
445:       if (rank == 0) {
446:         if (r) {
447:           PetscMPIInt nn;
448:           PetscCallMPI(MPI_Recv(array, nnodes * bs, MPIU_SCALAR, r, tag, comm, &status));
449:           PetscCallMPI(MPI_Get_count(&status, MPIU_SCALAR, &nn));
450:           PetscCheck(nn == nnodes * bs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Array size mismatch receiving from rank %" PetscInt_FMT, r);
451:         } else PetscCall(PetscArraycpy(array, x, nnodes * bs));
452:         /* If any fields are named, add scalar fields. Otherwise, add a vector field */
453:         PetscCall(DMDAGetFieldsNamed(daCurr, &fieldsnamed));
454:         if (fieldsnamed) {
455:           for (f = 0; f < bs; f++) {
456:             /* Extract and transpose the f'th field */
457:             for (k = 0; k < zm; k++) {
458:               for (j = 0; j < ym; j++) {
459:                 for (i = 0; i < xm; i++) {
460:                   PetscInt Iloc = i + xm * (j + ym * k);
461:                   array2[Iloc]  = array[Iloc * bs + f];
462:                 }
463:               }
464:             }
465:             PetscCall(PetscViewerVTKFWrite(viewer, fp, array2, nnodes, MPIU_SCALAR));
466:           }
467:         }
468:         PetscCall(PetscViewerVTKFWrite(viewer, fp, array, nnodes * bs, MPIU_SCALAR));

470:       } else if (r == rank) {
471:         PetscCallMPI(MPI_Send((void *)x, nnodes * bs, MPIU_SCALAR, 0, tag, comm));
472:       }
473:       PetscCall(VecRestoreArrayRead(X, &x));
474:     }
475:   }
476:   PetscCall(PetscFree2(array, array2));
477:   PetscCall(PetscFree(grloc));

479:   PetscCall(PetscFPrintf(comm, fp, "\n </AppendedData>\n"));
480:   PetscCall(PetscFPrintf(comm, fp, "</VTKFile>\n"));
481:   PetscCall(PetscFClose(comm, fp));
482:   PetscFunctionReturn(PETSC_SUCCESS);
483: }

485: /*@C
486:   DMDAVTKWriteAll - Write a file containing all the fields that have been provided to the viewer

488:   Collective

490:   Input Parameters:
491: + odm    - `DMDA` specifying the grid layout, passed as a `PetscObject`
492: - viewer - viewer of type `PETSCVIEWERVTK`

494:   Level: developer

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

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

504: .seealso: [](sec_struct), `DMDA`, `DM`, `PETSCVIEWERVTK`, `DMDASetFieldName()`
505: @*/
506: PetscErrorCode DMDAVTKWriteAll(PetscObject odm, PetscViewer viewer)
507: {
508:   DM        dm = (DM)odm;
509:   PetscBool isvtk;

511:   PetscFunctionBegin;
514:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk));
515:   PetscCheck(isvtk, PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_INCOMP, "Cannot use viewer type %s", ((PetscObject)viewer)->type_name);
516:   switch (viewer->format) {
517:   case PETSC_VIEWER_VTK_VTS:
518:     PetscCall(DMDAVTKWriteAll_VTS(dm, viewer));
519:     break;
520:   case PETSC_VIEWER_VTK_VTR:
521:     PetscCall(DMDAVTKWriteAll_VTR(dm, viewer));
522:     break;
523:   default:
524:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
525:   }
526:   PetscFunctionReturn(PETSC_SUCCESS);
527: }