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: const char *precision = PetscDefined(USE_REAL_SINGLE) ? "Float32" : (PetscDefined(USE_REAL_DOUBLE) ? "Float64" : "UnknownPrecision");
33: MPI_Comm comm;
34: Vec Coords;
35: PetscViewer_VTK *vtk = (PetscViewer_VTK *)viewer->data;
36: PetscViewerVTKObjectLink link;
37: FILE *fp;
38: PetscMPIInt rank, size, tag;
39: DMDALocalInfo info;
40: PetscInt dim, mx, my, mz, cdim, bs, maxnnodes, maxbs, i, j, k;
41: PetscInt64 boffset;
42: PetscInt rloc[6], (*grloc)[6] = NULL;
43: PetscScalar *array, *array2;
45: PetscFunctionBegin;
46: PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
47: PetscCheck(!PetscDefined(USE_COMPLEX), PETSC_COMM_SELF, PETSC_ERR_SUP, "Complex values not supported");
48: PetscCallMPI(MPI_Comm_size(comm, &size));
49: PetscCallMPI(MPI_Comm_rank(comm, &rank));
50: PetscCall(DMDAGetInfo(da, &dim, &mx, &my, &mz, NULL, NULL, NULL, &bs, NULL, NULL, NULL, NULL, NULL));
51: PetscCall(DMDAGetLocalInfo(da, &info));
52: PetscCall(DMGetCoordinates(da, &Coords));
53: if (Coords) {
54: PetscInt csize;
55: PetscCall(VecGetSize(Coords, &csize));
56: PetscCheck(csize % (mx * my * mz) == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Coordinate vector size mismatch");
57: cdim = csize / (mx * my * mz);
58: PetscCheck(cdim >= dim && cdim <= 3, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Coordinate vector size mismatch");
59: } else {
60: cdim = dim;
61: }
63: PetscCall(PetscFOpen(comm, vtk->filename, "wb", &fp));
64: PetscCall(PetscFPrintf(comm, fp, "<?xml version=\"1.0\"?>\n"));
65: PetscCall(PetscFPrintf(comm, fp, "<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"%s\" header_type=\"UInt64\">\n", byte_order));
66: 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));
68: if (rank == 0) PetscCall(PetscMalloc1(size * 6, &grloc));
69: rloc[0] = info.xs;
70: rloc[1] = info.xm;
71: rloc[2] = info.ys;
72: rloc[3] = info.ym;
73: rloc[4] = info.zs;
74: rloc[5] = info.zm;
75: PetscCallMPI(MPI_Gather(rloc, 6, MPIU_INT, rank == 0 ? grloc[0] : NULL, 6, MPIU_INT, 0, comm));
77: /* Write XML header */
78: maxnnodes = 0; /* Used for the temporary array size on rank 0 */
79: maxbs = 0; /* Used for the temporary array size on rank 0 */
80: boffset = 0; /* Offset into binary file */
81: for (PetscMPIInt r = 0; r < size; r++) {
82: PetscInt xs = -1, xm = -1, ys = -1, ym = -1, zs = -1, zm = -1, nnodes = 0;
84: if (rank == 0) {
85: xs = grloc[r][0];
86: xm = grloc[r][1];
87: ys = grloc[r][2];
88: ym = grloc[r][3];
89: zs = grloc[r][4];
90: zm = grloc[r][5];
91: nnodes = xm * ym * zm;
92: }
93: maxnnodes = PetscMax(maxnnodes, nnodes);
94: 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));
95: PetscCall(PetscFPrintf(comm, fp, " <Points>\n"));
96: PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"Position\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, boffset));
97: boffset += 3 * nnodes * sizeof(PetscScalar) + sizeof(PetscInt64);
98: PetscCall(PetscFPrintf(comm, fp, " </Points>\n"));
100: PetscCall(PetscFPrintf(comm, fp, " <PointData Scalars=\"ScalarPointData\">\n"));
101: for (link = vtk->link; link; link = link->next) {
102: Vec X = (Vec)link->vec;
103: PetscInt bs;
104: DM daCurr;
105: PetscBool fieldsnamed;
106: const char *vecname = "Unnamed";
108: PetscCall(VecGetDM(X, &daCurr));
109: PetscCall(DMDAGetInfo(daCurr, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &bs, NULL, NULL, NULL, NULL, NULL));
110: maxbs = PetscMax(maxbs, bs);
112: if (((PetscObject)X)->name || link != vtk->link) PetscCall(PetscObjectGetName((PetscObject)X, &vecname));
114: /* If any fields are named, add scalar fields. Otherwise, add a vector field */
115: PetscCall(DMDAGetFieldsNamed(daCurr, &fieldsnamed));
116: if (fieldsnamed) {
117: for (PetscInt f = 0; f < bs; f++) {
118: char buf[256];
119: const char *fieldname;
120: PetscCall(DMDAGetFieldName(daCurr, f, &fieldname));
121: if (!fieldname) {
122: PetscCall(PetscSNPrintf(buf, sizeof(buf), "%" PetscInt_FMT, f));
123: fieldname = buf;
124: }
125: PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s.%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset));
126: boffset += nnodes * sizeof(PetscScalar) + sizeof(PetscInt64);
127: }
128: } else {
129: PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s\" NumberOfComponents=\"%" PetscInt_FMT "\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, bs, boffset));
130: boffset += bs * nnodes * sizeof(PetscScalar) + sizeof(PetscInt64);
131: }
132: }
133: PetscCall(PetscFPrintf(comm, fp, " </PointData>\n"));
134: PetscCall(PetscFPrintf(comm, fp, " </Piece>\n"));
135: }
136: PetscCall(PetscFPrintf(comm, fp, " </StructuredGrid>\n"));
137: PetscCall(PetscFPrintf(comm, fp, " <AppendedData encoding=\"raw\">\n"));
138: PetscCall(PetscFPrintf(comm, fp, "_"));
140: /* Now write the arrays. */
141: tag = ((PetscObject)viewer)->tag;
142: PetscCall(PetscMalloc2(maxnnodes * PetscMax(3, maxbs), &array, maxnnodes * PetscMax(3, maxbs), &array2));
143: for (PetscMPIInt r = 0; r < size; r++) {
144: MPI_Status status;
145: PetscInt xs = -1, xm = -1, ys = -1, ym = -1, zs = -1, zm = -1, nnodes = 0;
147: if (rank == 0) {
148: xs = grloc[r][0];
149: xm = grloc[r][1];
150: ys = grloc[r][2];
151: ym = grloc[r][3];
152: zs = grloc[r][4];
153: zm = grloc[r][5];
154: nnodes = xm * ym * zm;
155: } else if (r == rank) {
156: nnodes = info.xm * info.ym * info.zm;
157: }
159: /* Write the coordinates */
160: if (Coords) {
161: const PetscScalar *coords;
163: PetscCall(VecGetArrayRead(Coords, &coords));
164: if (rank == 0) {
165: if (r) {
166: PetscMPIInt nn;
168: PetscCallMPI(MPIU_Recv(array, nnodes * cdim, MPIU_SCALAR, r, tag, comm, &status));
169: PetscCallMPI(MPI_Get_count(&status, MPIU_SCALAR, &nn));
170: PetscCheck(nn == nnodes * cdim, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Array size mismatch");
171: } else PetscCall(PetscArraycpy(array, coords, nnodes * cdim));
172: /* Transpose coordinates to VTK (C-style) ordering */
173: for (k = 0; k < zm; k++) {
174: for (j = 0; j < ym; j++) {
175: for (i = 0; i < xm; i++) {
176: PetscInt Iloc = i + xm * (j + ym * k);
177: array2[Iloc * 3 + 0] = array[Iloc * cdim + 0];
178: array2[Iloc * 3 + 1] = cdim > 1 ? array[Iloc * cdim + 1] : 0.0;
179: array2[Iloc * 3 + 2] = cdim > 2 ? array[Iloc * cdim + 2] : 0.0;
180: }
181: }
182: }
183: } else if (r == rank) {
184: PetscCallMPI(MPIU_Send((void *)coords, nnodes * cdim, MPIU_SCALAR, 0, tag, comm));
185: }
186: PetscCall(VecRestoreArrayRead(Coords, &coords));
187: } else { /* Fabricate some coordinates using grid index */
188: for (k = 0; k < zm; k++) {
189: for (j = 0; j < ym; j++) {
190: for (i = 0; i < xm; i++) {
191: PetscInt Iloc = i + xm * (j + ym * k);
192: array2[Iloc * 3 + 0] = xs + i;
193: array2[Iloc * 3 + 1] = ys + j;
194: array2[Iloc * 3 + 2] = zs + k;
195: }
196: }
197: }
198: }
199: PetscCall(PetscViewerVTKFWrite(viewer, fp, array2, nnodes * 3, MPIU_SCALAR));
201: /* Write each of the objects queued up for this file */
202: for (link = vtk->link; link; link = link->next) {
203: Vec X = (Vec)link->vec;
204: const PetscScalar *x;
205: PetscInt bs;
206: DM daCurr;
207: PetscBool fieldsnamed;
208: PetscCall(VecGetDM(X, &daCurr));
209: PetscCall(DMDAGetInfo(daCurr, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &bs, NULL, NULL, NULL, NULL, NULL));
210: PetscCall(VecGetArrayRead(X, &x));
211: if (rank == 0) {
212: if (r) {
213: PetscMPIInt nn;
215: PetscCallMPI(MPIU_Recv(array, nnodes * bs, MPIU_SCALAR, r, tag, comm, &status));
216: PetscCallMPI(MPI_Get_count(&status, MPIU_SCALAR, &nn));
217: PetscCheck(nn == nnodes * bs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Array size mismatch receiving from rank %d", r);
218: } else PetscCall(PetscArraycpy(array, x, nnodes * bs));
220: /* If any fields are named, add scalar fields. Otherwise, add a vector field */
221: PetscCall(DMDAGetFieldsNamed(daCurr, &fieldsnamed));
222: if (fieldsnamed) {
223: for (PetscInt f = 0; f < bs; f++) {
224: /* Extract and transpose the f'th field */
225: for (k = 0; k < zm; k++) {
226: for (j = 0; j < ym; j++) {
227: for (i = 0; i < xm; i++) {
228: PetscInt Iloc = i + xm * (j + ym * k);
229: array2[Iloc] = array[Iloc * bs + f];
230: }
231: }
232: }
233: PetscCall(PetscViewerVTKFWrite(viewer, fp, array2, nnodes, MPIU_SCALAR));
234: }
235: } else PetscCall(PetscViewerVTKFWrite(viewer, fp, array, bs * nnodes, MPIU_SCALAR));
236: } else if (r == rank) PetscCallMPI(MPIU_Send((void *)x, nnodes * bs, MPIU_SCALAR, 0, tag, comm));
237: PetscCall(VecRestoreArrayRead(X, &x));
238: }
239: }
240: PetscCall(PetscFree2(array, array2));
241: PetscCall(PetscFree(grloc));
243: PetscCall(PetscFPrintf(comm, fp, "\n </AppendedData>\n"));
244: PetscCall(PetscFPrintf(comm, fp, "</VTKFile>\n"));
245: PetscCall(PetscFClose(comm, fp));
246: PetscFunctionReturn(PETSC_SUCCESS);
247: }
249: static PetscErrorCode DMDAVTKWriteAll_VTR(DM da, PetscViewer viewer)
250: {
251: const char *byte_order = PetscBinaryBigEndian() ? "BigEndian" : "LittleEndian";
252: const char *precision = PetscDefined(USE_REAL_SINGLE) ? "Float32" : (PetscDefined(USE_REAL_DOUBLE) ? "Float64" : "UnknownPrecision");
253: MPI_Comm comm;
254: PetscViewer_VTK *vtk = (PetscViewer_VTK *)viewer->data;
255: PetscViewerVTKObjectLink link;
256: FILE *fp;
257: PetscMPIInt rank, size, tag;
258: DMDALocalInfo info;
259: PetscInt dim, mx, my, mz, maxnnodes, maxbs, i, j, k;
260: PetscInt64 boffset;
261: PetscInt rloc[6], (*grloc)[6] = NULL;
262: PetscScalar *array, *array2;
264: PetscFunctionBegin;
265: PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
266: PetscCheck(!PetscDefined(USE_COMPLEX), PETSC_COMM_SELF, PETSC_ERR_SUP, "Complex values not supported");
267: PetscCallMPI(MPI_Comm_size(comm, &size));
268: PetscCallMPI(MPI_Comm_rank(comm, &rank));
269: PetscCall(DMDAGetInfo(da, &dim, &mx, &my, &mz, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
270: PetscCall(DMDAGetLocalInfo(da, &info));
271: PetscCall(PetscFOpen(comm, vtk->filename, "wb", &fp));
272: PetscCall(PetscFPrintf(comm, fp, "<?xml version=\"1.0\"?>\n"));
273: PetscCall(PetscFPrintf(comm, fp, "<VTKFile type=\"RectilinearGrid\" version=\"0.1\" byte_order=\"%s\" header_type=\"UInt64\">\n", byte_order));
274: 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));
276: if (rank == 0) PetscCall(PetscMalloc1(size * 6, &grloc));
277: rloc[0] = info.xs;
278: rloc[1] = info.xm;
279: rloc[2] = info.ys;
280: rloc[3] = info.ym;
281: rloc[4] = info.zs;
282: rloc[5] = info.zm;
283: PetscCallMPI(MPI_Gather(rloc, 6, MPIU_INT, rank == 0 ? grloc[0] : NULL, 6, MPIU_INT, 0, comm));
285: /* Write XML header */
286: maxnnodes = 0; /* Used for the temporary array size on rank 0 */
287: maxbs = 0; /* Used for the temporary array size on rank 0 */
288: boffset = 0; /* Offset into binary file */
289: for (PetscMPIInt r = 0; r < size; r++) {
290: PetscInt xs = -1, xm = -1, ys = -1, ym = -1, zs = -1, zm = -1, nnodes = 0;
292: if (rank == 0) {
293: xs = grloc[r][0];
294: xm = grloc[r][1];
295: ys = grloc[r][2];
296: ym = grloc[r][3];
297: zs = grloc[r][4];
298: zm = grloc[r][5];
299: nnodes = xm * ym * zm;
300: }
301: maxnnodes = PetscMax(maxnnodes, nnodes);
302: 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));
303: PetscCall(PetscFPrintf(comm, fp, " <Coordinates>\n"));
304: PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"Xcoord\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, boffset));
305: boffset += xm * sizeof(PetscScalar) + sizeof(PetscInt64);
306: PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"Ycoord\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, boffset));
307: boffset += ym * sizeof(PetscScalar) + sizeof(PetscInt64);
308: PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"Zcoord\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, boffset));
309: boffset += zm * sizeof(PetscScalar) + sizeof(PetscInt64);
310: PetscCall(PetscFPrintf(comm, fp, " </Coordinates>\n"));
311: PetscCall(PetscFPrintf(comm, fp, " <PointData Scalars=\"ScalarPointData\">\n"));
312: for (link = vtk->link; link; link = link->next) {
313: Vec X = (Vec)link->vec;
314: PetscInt bs;
315: DM daCurr;
316: PetscBool fieldsnamed;
317: const char *vecname = "Unnamed";
319: PetscCall(VecGetDM(X, &daCurr));
320: PetscCall(DMDAGetInfo(daCurr, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &bs, NULL, NULL, NULL, NULL, NULL));
321: maxbs = PetscMax(maxbs, bs);
322: if (((PetscObject)X)->name || link != vtk->link) PetscCall(PetscObjectGetName((PetscObject)X, &vecname));
324: /* If any fields are named, add scalar fields. Otherwise, add a vector field */
325: PetscCall(DMDAGetFieldsNamed(daCurr, &fieldsnamed));
326: if (fieldsnamed) {
327: for (PetscInt f = 0; f < bs; f++) {
328: char buf[256];
329: const char *fieldname;
330: PetscCall(DMDAGetFieldName(daCurr, f, &fieldname));
331: if (!fieldname) {
332: PetscCall(PetscSNPrintf(buf, sizeof(buf), "%" PetscInt_FMT, f));
333: fieldname = buf;
334: }
335: PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s.%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset));
336: boffset += nnodes * sizeof(PetscScalar) + sizeof(PetscInt64);
337: }
338: } else {
339: PetscCall(PetscFPrintf(comm, fp, " <DataArray type=\"%s\" Name=\"%s\" NumberOfComponents=\"%" PetscInt_FMT "\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, bs, boffset));
340: boffset += bs * nnodes * sizeof(PetscScalar) + sizeof(PetscInt64);
341: }
342: }
343: PetscCall(PetscFPrintf(comm, fp, " </PointData>\n"));
344: PetscCall(PetscFPrintf(comm, fp, " </Piece>\n"));
345: }
346: PetscCall(PetscFPrintf(comm, fp, " </RectilinearGrid>\n"));
347: PetscCall(PetscFPrintf(comm, fp, " <AppendedData encoding=\"raw\">\n"));
348: PetscCall(PetscFPrintf(comm, fp, "_"));
350: /* Now write the arrays. */
351: tag = ((PetscObject)viewer)->tag;
352: PetscCall(PetscMalloc2(PetscMax(maxnnodes * maxbs, info.xm + info.ym + info.zm), &array, PetscMax(maxnnodes * maxbs, info.xm + info.ym + info.zm), &array2));
353: for (PetscMPIInt r = 0; r < size; r++) {
354: MPI_Status status;
355: PetscInt xs = -1, xm = -1, ys = -1, ym = -1, zs = -1, zm = -1, nnodes = 0;
357: if (rank == 0) {
358: xs = grloc[r][0];
359: xm = grloc[r][1];
360: ys = grloc[r][2];
361: ym = grloc[r][3];
362: zs = grloc[r][4];
363: zm = grloc[r][5];
364: nnodes = xm * ym * zm;
365: } else if (r == rank) {
366: nnodes = info.xm * info.ym * info.zm;
367: }
368: { /* Write the coordinates */
369: Vec Coords;
371: PetscCall(DMGetCoordinates(da, &Coords));
372: if (Coords) {
373: const PetscScalar *coords;
374: PetscCall(VecGetArrayRead(Coords, &coords));
375: if (rank == 0) {
376: if (r) {
377: PetscMPIInt nn;
379: PetscCallMPI(MPIU_Recv(array, xm + ym + zm, MPIU_SCALAR, r, tag, comm, &status));
380: PetscCallMPI(MPI_Get_count(&status, MPIU_SCALAR, &nn));
381: PetscCheck(nn == xm + ym + zm, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Array size mismatch");
382: } else {
383: /* x coordinates */
384: for (j = 0, k = 0, i = 0; i < xm; i++) {
385: PetscInt Iloc = i + xm * (j + ym * k);
387: array[i] = coords[Iloc * dim + 0];
388: }
389: /* y coordinates */
390: for (i = 0, k = 0, j = 0; j < ym; j++) {
391: PetscInt Iloc = i + xm * (j + ym * k);
393: array[j + xm] = dim > 1 ? coords[Iloc * dim + 1] : 0;
394: }
395: /* z coordinates */
396: for (i = 0, j = 0, k = 0; k < zm; k++) {
397: PetscInt Iloc = i + xm * (j + ym * k);
399: array[k + xm + ym] = dim > 2 ? coords[Iloc * dim + 2] : 0;
400: }
401: }
402: } else if (r == rank) {
403: xm = info.xm;
404: ym = info.ym;
405: zm = info.zm;
406: for (j = 0, k = 0, i = 0; i < xm; i++) {
407: 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);
414: array2[j + xm] = dim > 1 ? coords[Iloc * dim + 1] : 0;
415: }
416: for (i = 0, j = 0, k = 0; k < zm; k++) {
417: PetscInt Iloc = i + xm * (j + ym * k);
419: array2[k + xm + ym] = dim > 2 ? coords[Iloc * dim + 2] : 0;
420: }
421: PetscCallMPI(MPIU_Send((void *)array2, xm + ym + zm, MPIU_SCALAR, 0, tag, comm));
422: }
423: PetscCall(VecRestoreArrayRead(Coords, &coords));
424: } else { /* Fabricate some coordinates using grid index */
425: for (i = 0; i < xm; i++) array[i] = xs + i;
426: for (j = 0; j < ym; j++) array[j + xm] = ys + j;
427: for (k = 0; k < zm; k++) array[k + xm + ym] = zs + k;
428: }
429: if (rank == 0) {
430: PetscCall(PetscViewerVTKFWrite(viewer, fp, &array[0], xm, MPIU_SCALAR));
431: PetscCall(PetscViewerVTKFWrite(viewer, fp, &array[xm], ym, MPIU_SCALAR));
432: PetscCall(PetscViewerVTKFWrite(viewer, fp, &array[xm + ym], zm, MPIU_SCALAR));
433: }
434: }
436: /* Write each of the objects queued up for this file */
437: for (link = vtk->link; link; link = link->next) {
438: Vec X = (Vec)link->vec;
439: const PetscScalar *x;
440: PetscInt bs;
441: DM daCurr;
442: PetscBool fieldsnamed;
444: PetscCall(VecGetDM(X, &daCurr));
445: PetscCall(DMDAGetInfo(daCurr, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &bs, NULL, NULL, NULL, NULL, NULL));
447: PetscCall(VecGetArrayRead(X, &x));
448: if (rank == 0) {
449: if (r) {
450: PetscMPIInt nn;
452: PetscCallMPI(MPIU_Recv(array, nnodes * bs, MPIU_SCALAR, r, tag, comm, &status));
453: PetscCallMPI(MPI_Get_count(&status, MPIU_SCALAR, &nn));
454: PetscCheck(nn == nnodes * bs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Array size mismatch receiving from rank %d", r);
455: } else PetscCall(PetscArraycpy(array, x, nnodes * bs));
456: /* If any fields are named, add scalar fields. Otherwise, add a vector field */
457: PetscCall(DMDAGetFieldsNamed(daCurr, &fieldsnamed));
458: if (fieldsnamed) {
459: for (PetscInt f = 0; f < bs; f++) {
460: /* Extract and transpose the f'th field */
461: for (k = 0; k < zm; k++) {
462: for (j = 0; j < ym; j++) {
463: for (i = 0; i < xm; i++) {
464: PetscInt Iloc = i + xm * (j + ym * k);
465: array2[Iloc] = array[Iloc * bs + f];
466: }
467: }
468: }
469: PetscCall(PetscViewerVTKFWrite(viewer, fp, array2, nnodes, MPIU_SCALAR));
470: }
471: } else PetscCall(PetscViewerVTKFWrite(viewer, fp, array, nnodes * bs, MPIU_SCALAR));
472: } else if (r == rank) {
473: PetscCallMPI(MPIU_Send((void *)x, nnodes * bs, MPIU_SCALAR, 0, tag, comm));
474: }
475: PetscCall(VecRestoreArrayRead(X, &x));
476: }
477: }
478: PetscCall(PetscFree2(array, array2));
479: PetscCall(PetscFree(grloc));
481: PetscCall(PetscFPrintf(comm, fp, "\n </AppendedData>\n"));
482: PetscCall(PetscFPrintf(comm, fp, "</VTKFile>\n"));
483: PetscCall(PetscFClose(comm, fp));
484: PetscFunctionReturn(PETSC_SUCCESS);
485: }
487: /*@
488: DMDAVTKWriteAll - Write a file containing all the fields that have been provided to the viewer
490: Collective
492: Input Parameters:
493: + odm - `DMDA` specifying the grid layout, passed as a `PetscObject`
494: - viewer - viewer of type `PETSCVIEWERVTK`
496: Level: developer
498: Notes:
499: This function is a callback used by the `PETSCVIEWERVTK` viewer to actually write the file.
500: The reason for this odd model is that the VTK file format does not provide any way to write one field at a time.
501: Instead, metadata for the entire file needs to be available up-front before you can start writing the file.
503: If any fields have been named (see e.g. `DMDASetFieldName()`), then individual scalar
504: fields are written. Otherwise, a single multi-dof (vector) field is written.
506: .seealso: [](sec_struct), `DMDA`, `DM`, `PETSCVIEWERVTK`, `DMDASetFieldName()`
507: @*/
508: PetscErrorCode DMDAVTKWriteAll(PetscObject odm, PetscViewer viewer)
509: {
510: DM dm = (DM)odm;
511: PetscBool isvtk;
513: PetscFunctionBegin;
516: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk));
517: PetscCheck(isvtk, PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_INCOMP, "Cannot use viewer type %s", ((PetscObject)viewer)->type_name);
518: switch (viewer->format) {
519: case PETSC_VIEWER_VTK_VTS:
520: PetscCall(DMDAVTKWriteAll_VTS(dm, viewer));
521: break;
522: case PETSC_VIEWER_VTK_VTR:
523: PetscCall(DMDAVTKWriteAll_VTR(dm, viewer));
524: break;
525: default:
526: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
527: }
528: PetscFunctionReturn(PETSC_SUCCESS);
529: }