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: PetscCall(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(MPI_Recv(&numCorners, 1, MPIU_INT, proc, tag, comm, &status));
156: PetscCallMPI(MPI_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(MPI_Send(&numSend, 1, MPIU_INT, 0, tag, comm));
195: PetscCallMPI(MPI_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(MPI_Recv(&numCells, 1, MPIU_INT, proc, tag, comm, &status));
211: PetscCallMPI(MPI_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(MPI_Send(&numCells, 1, MPIU_INT, 0, tag, comm));
219: PetscCallMPI(MPI_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: PetscCall(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(MPI_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(MPI_Send(&k, 1, MPIU_INT, 0, tag, comm));
400: PetscCallMPI(MPI_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: PetscCall(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, §ion)); /* 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 *)§ion));
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, §ion));
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, §ion));
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 *)§ion));
597: if (!section) {
598: DM dmX;
600: PetscCall(VecGetDM(X, &dmX));
601: if (dmX) PetscCall(DMGetLocalSection(dmX, §ion));
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, §ion));
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: /*@C
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: }